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ABSTRACT 



The search for a universal solution of the equations of motion for a satellite or- 
biting an oblate planet is a subject that has merited great interest because of its 
theoretical implications and practical applications. The discovery of such a solu- 
tion should motivate a reassessment of both the theories that exhibit singularities 
and the physical effects implied by singularities. The practical importance of such 
a solution is the efficiency of simple analytic formulas in predicting simultaneously 
the paths of large numbers of satellites in a multitude of orbits. Here, a complete 
first order solution to the problem of a satellite, perturbed only by the oblateness 
of the Earth, is displayed. The orbit is free of singularities for all parameters and 
is valid for 1000 revolutions with a relative error of the order « 10“®. 
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NOTATION 



a semi major axis of the initial instantaneous ellipse 
c cos to 

e eccentricity of the initial instaintaneous ellipse (0 < e << l) 
h constant value of ^ along am orbit, approximately equal to the 

initial magnitute of angular momentum 
i inclination of reference plane 
t'o initial value of i 

J 2 oblateness coefficient of the planet (coefficient of the second harmonic in 
the expansion of the gravitational potential) 




p I GM where G is the gravitational constant and M is the mass of the 

planet 

R equatorial radius of the planet 

r radial distance from the center of the planet to the satellite 
s sin I'o 
t time 

to initial value of t 
u 2 

r 

V gravitational potential 

6 latitude of the satellite 

<j> longitude of the satellite 

6 angle from line NON' to the satellite measured in the reference plane 
(Fig. 5.1) 

6q initial value of 6 
n longitude of line NON' (Fig. 5.1) 
flo initial value of 0 

U-’ argument of perigee of the initial instantaneous ellipse 
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I. INTRODUCTION 



A characteristic feature of practical satellite orbit prediction is that the engi- 
neer may deal with numerous satellites in a great variety of alternative orbits. Un- 
der these and many other such circumstances analytic relations which can quickly 
approximate an orbit may be far superior to large numerical models. While many 
analytic methods have been developed for the artificial satellite age, most are not 
used in practical orbit prediction because they violate one or more of the following 
principles: 

• The method should provide a solution that is significantly more accurate 
than the two-body solution. 

• The real physical effects of the orbit should be easily distinguishable in the 
solution. 

• The solution should be universal; it should be valid for all orbital parameters. 

The motivation for this research was the desire to develop a method for satellite 
prediction that would embody these characteristics. 

In this analysis, a solution to the equations of motion of a satellite around an 
oblate planet is found by use of a variation of the perturbation technique known 
as the Method of Strained Coordinates. The orbit is valid for 1000 revolutions 
with a relative error of 10“®. The solution which is valid for all eccentricities and 
for all inclinations, was obtained by extensive use of the symbolic manipulation 
program MACSYMA. 

The analysis begins with a background discussion of some of the competing 
satellite orbit theories. There is then a development of the equations of motion 
beginning with a deri\'ation of the two-body solution. The various forces which act 
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to disturb the two-body orbit are highlighted; a more thorough discussion is given 
for the effects of oblateness. There is a complete treatment of the perturbation 
technique as the equations of motion are solved in detail. The complete first order 
solution is displayed as a function of coordinates and as a function of the orbital 
elements. In addition, an independent analysis of the polar and equatorial orbits 
is performed to serve as a check of the general solution. 
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II. BACKGROUND 



A. INTRODUCTION 

The theory of flight of artificial satellites is closely related to claissical celes- 
tial mechajiics, one of the oldest and most highly developed branches of science. 
The equations which describe the motion of artificial satellite are in principle 
identical to the equations of motion of natural celestial bodies. It is not surprising 
then that the results originally derived in classical celestial mechanics have been 
freely used to explain the motion of artificial celestial bodies. 

The foundations of classical celestial mechanics were established in the eigh- 
teenth century when Clairaut, d’Alembert, Laplace, Lagrange, and Euler intro- 
duced theories and analytical methods to explain the large deviation of the Moon 
from an elliptic orbit due to solar attraction. These theories all supplemented 
or complemented the pioneer work that had been done by Newton. Newton had 
correctly indicated the Moon’s variation in eccentricity and inclination and the re- 
gression of the nodes; however, his published theory accounted for only half of the 
motion of perigee. (A century later an unpublished work was found to contain the 
full explanation.) Clairaut in 1749 was on the verge of substituting a new law of 
gravitation for the Newtonian law when he found that second order perturbations 
removed the discrepancy in the motion of perigee. 

In that same century, Euler investigated and developed the perturbative 
function and began the development of the method of variation of parameters, 
which Wtis later extended by Lagrange. 
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Lagrange, Laplace, and Poisson all advanced the discipline through their 
investigation of the stability of the solar system. Later, in the nineteenth cen- 
tury, the use of Hamiltonian mechanics was used to great advantage by Delaunay 
whose work influenced the innovators of artificial satellite orbit prediction, namely 
Brouwer, Kozai, and GMfinkel. 

The introduction of the modern computer and the launch of the first artifi- 
cial satellites had a profound effect on the science of celestial mechanics. Classical 
celestial mechanics had been essentially a contemplative science with the principle 
aim being to study the laws of motion of existing heavenly bodies. In contrast, the 
science of the flight of artificial satellites is an active engineering science concerned 
with determining or predicting relatively short term orbits, and in many cases con- 
trolling the satellite’s motion through on-board propulsive devices. It was in the 
exciting climate, following the successful launchings of the first artificial satellites, 
that most of the new methods of satellite orbit prediction were developed. 

B. THE USE OF HAMILTONIAN MECHANICS IN SATELLITE 
ORBIT PREDICTION 

There are certain schools of thought in dynamic astronomy and theoretical 
physics that support the loyal use of Hamiltonian mechanics [Ref. 1, p. 228], 
and many of the new methods in general perturbations take advantage of the 
elegant formalism offered by the Hamiltonian method. The Hamiltonian method 
is referred to here as the formal process of writing the equations of motion for a 
satellite in the canonical form: 



dqr dH 
dt dpr 



^ _ -dH 
dt dqr 



(r = l,2...3n) 



where q^ = generalized coordinates 



( 2 . 1 ) 
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Pr = generalized momenta 
H - r 



{H is the Hamiltonian and L is the Lagrangian, L = T — V , kinetic energy - 
potential energy). 

The solution to (2.1) may be written down if a function S can be found, where S 
is any complete solution of the Hamilton- Jau:obi equation 




It should be noted that (2.2) is tractable only if the variables gr and t are separable 
within 5 and H. 

Although this present analysis does not use Hamiltonian mechanics to solve 
the equations of motion, a summary of its use is merited since the advances in this 
area have been an invaluable contribution to general perturbation theory. 

In celestial mechanics, one may formulate a Hamiltonian that represents the 
gravitational attraction of the central force and add to it terms which are “per- 
turbing Hamiltonians.” Various sets of canonic variables may be chosen with the 
goal of expressing a zero-order Hamiltonian in a simple form and the higher order 
effects in an iterative fashion. Each term then may be dealt with through a succes- 
sion of canonic transformations. Delauney introduced a systematic procedure for 
isolating parts of the Hamiltonian and then generating a suitable transformation 
in successive steps. A particular feature of his approach is that a periodic term in 
the Hamiltonian may be eliminated with each canonical transformation [Ref. 2]. 
While DeLauney used a procedure that eliminated one term at a time, von Ziepel 
devised a technique that eliminates one angular variable with each transformation. 
This method reduces the number of degrees of freedom while at the same time 
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impjLTting to the transformed Hamiltonian a symmetry of the unperturbed system 
[Ref. l]. The eliminated variables are referred to as ignorable coordinates since 
they do not participate in the solution of the transformed equations of motion but 
can be recovered after the solution has been obtained. 

Using von Ziepel’s method, Brouwer devised one of the most notable general 
perturbation theories [Ref. 3]. Prior to Brouwer’s method, all previous work in 
artificial satellite theory had written the Hamiltonian as a Fourier series in the 
mean anomaly with coefficients that were infinite series in powers of the eccen- 
tricity. Brouwer used an elliptic approximation for the potential and obtained a 
complete first order theory with some second order development using canonical 
transformations. The essence of Brouwer’s method was to write equations (2.1) as 



^ ^ _ dH 

dt dlj dt dLj 



where 



Li = VG M^ h=M 

Li — Li\/l — li = u! 

Lz = Li cos t I 3 = Cl 



are the Delauney variables where 



a - semimajor axis 
i - inclination 
w - argument of perigee 



€ - eccentricity 
M - mean anomaly 
fl - longitude of 

ascending node 



(2.3) 



Referring back to the Hamilton- Jacobi equation (2.2), the function S is used as 
a generating function to find a new Hamiltonian that leads to simplified canonic 
equations. By choosing S correctly, Brouwer was able to find a canonical transfor- 
mation from the Delauney variables to a set of double primed Delauney variables 
(L", I'-) such that (2.3) have the form 
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^ _ 0 

dt ~ dt ~ dLj 

where H’ is the Hamiltonian expressed in terms of double primed orbital elements 

by 

V[ = \fGM^ I'l = M" 

L'j' = Li'vT^^ l'i = w" 

L'l = L'icosi" l'^ = n" 

The double primed orbital elements are related to the unprimed elements by 

a = a" + ea e = t" + ce 
i = t" + et M = M" + cM 

w = w'* + tw n = n" + cn 

where the quantities €.a,€e,ei,eM,€w,eQ, are periodic functions of the double 
primed orbital elements. 

Equations (2.4) are solved for the double-primed Delauney variables which 
can then be expressed in terms of the original variables. 

The results initially obtained by Brouwer were not valid at the critical incli- 
nation of 63°.4, and they were questionable when either the inclination or eccen- 
tricity were near zero [Ref. 4). O. K. Smith devised a method for dealing with the 
problem of zero inclination and eccentricity [Ref. 5], but Brouwer subsequently 
challenged the validity of his method. Later, Lyddane [Ref. 6] was able to suc- 
cessfully remove the restrictions on small values of eccentricity and inclination by 
reformulating Brouwer’s work in terms of an alternate set of variables. 

Brouwer used the central force term as the first approximation for the po- 
tential since there is no exact solution to the equations of motion of a satellite 
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under the influence of the more complex potential described by an oblate planet 
having axial symmetry. Other authors have attempted to introduce a new poten- 
tial which approximates the Earth’s potential better than the central force term 
alone and also leads to an exact solution. Most notable hzis been the work of 
Sterne [Ref. 7) and Garfinkel [Ref. 8]. Sterne’s potential function accounts for 
most of the standard potential through the second harmonic and it leads to a 
solution of canonical constants that are free of first order secular perturbations. 
The remaining effects of the earth’s oblateness and other forces are allowed for 
in the perturbing Hamiltonian which causes the six canonical constants of the 
unperturbed solution to undergo variations with time. Sterne provided the inspi- 
ration for Garfinkel’s method which is essentially the same as Sterne’s but more 
developed. Garfinkel included the second and fourth harmonics and arrived at 
a solution that is reducible to quadratures. Garfinkel’s original solution was not 
valid at the critical inclination, but in a later paper he removed the singularity 
through a variation on his perturbation theory [Ref. 9). 

While Garfinkel did his analysis in spherical coordinates, Vinti [Ref. 10], 
derived a potential expressible in oblate spheroidal coordinates. In his original 
analysis, Vinti introduced a potential function and associated coordinate system 
that would lead to the separability of the Hamilton-Jacobi equations. In a subse- 
quent work [Ref. 11], Vinti showed that the equations result in a solution reducible 
to quadratures. The Vinti potential is an exact expression of the Earth’s potential 
through the second harmonic. The theory provides perturbed coordinates, not 
perturbed elements. In his second order analysis [Ref. 12], Kozai criticized this 
characteristic, and he chose not to use the Vinti potential since it would require 
changing the definitions of the conventional orbital elements. 
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Morrison [Ref. 13] showed that the von Zeipel method is a particular case of 
the method of averaging, and Liu [Ref. 14] used the latter technique to study the 
combined effects of air drag and planet oblateness on a satellite orbit. The method 
of averaging, unlike von Ziepel’s method, does not require that transformations 
be canonical. The method of averaging has been used extensively in recent years 
by Lorell and Liu [Ref. 15], McClain [Ref. 16], and Hoots [Ref. 17]. However, the 
validity of the method has been challenged; most notably Taff [Ref. 18] doubts 
its rigorous foundation. Arnold [Ref. 19] notes that the principle of averaging is a 
vaguely formulated and rigorously untrue assertion, but he adds that sometimes 
such assertions are fruitful mathematical sources. It should be also noted that 
the solutions obtained in Liu’s analysis [Ref. 14] are not valid for near circular 
equatorial orbits nor at the critical inclination, while those obtained by Hoots are 
valid only for small eccentricities (0 < c < .1), and they are invalid at inclinations 
near 0° or the critical. 

Hamiltonian mechanics has provided a rich source of literature in orbit the- 
ory; however its practical applicability has been questioned in some textbooks. 
Roy [Ref. 20] briefly outlines the use of the canonic equations, while Taff [Ref. 
18p. 322) states that he does not see any additional practical applications provided. 
Baker [Ref. 21] chooses not to represent the subject. A general criticism is that 
the process of generating suitable transformations in the perturbation procedure 
tends to make the coordinates and the momenta less distinguishable on physical 
grounds and more difficult to relate to the set of natural coordinates which were 
used to write down the initial set of differential equations. While the ultimate 
form of the governing equations may be simple to solve, there remains the tedious 
task of obtaining explicit results in terms of physically meaningful coordinates or 
elements. 
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C. KOZAI’S METHOD 



Prior to his work cited above, Kozai [Ref. 22] developed a method for find- 
ing the perturbations on the orbital elements of a satellite considering only the 
oblateness of the Earth. Kozai developed a disturbing function based on the 
Earth’s departure from a sphere, and he used a version of Lagrange’s planetary 
equations to formulate the solution. Kozai \ised the standard form for the Earth’s 
potential and included the harmonics and J 4 . Despite the use of the higher 

harmonics, the theory is first order. Kozai expressed short-period terms in Jj, the 
secular in Jj, J4, and J|, and the long period terms in Jj, and The analytic 
expressions are developed using the standard orbital elements. 

Kozai’s work is cited here because, due to its simplicity, the method has be- 
come very popular in many textbooks and handbooks on orbital mechanics. [Refs. 
20, 23, 24]. However, Taff cites Kozai’s method as an example of misapplication 
of perturbation theory [Ref. 18p. 332]. Taff challenges the assumptions made by 
Kozai in his analysis, and he points out that the method is invalid at the critical 
inclination. 

As was stated in Chapter I, a motivation for this current analysis was the 
purpose of finding a perturbation method that would lead to universal solution. All 
of the methods discussed thus far have particular problems at certain inclinations 
or eccentricities. Some of the problem cases have been resolved by unique efforts 
(These cases are discussed in Appendix C.); however, one should question the 
underlying validity of any perturbation method that produces singularities. There 
has been no satisfactory way found for avoiding the critical inclination singularity 
in Kozai’s method. 
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D. THE DIRECT METHODS 



R. E. Roberson [Ref. 25] devised an approach for finding the qualitative 
zmd approximate quantitative results concerning the behavior of a set of orbital 
elements in the gravitational field of an oblate plzmet. The motion of a neair 
satellite around the planet is simply described to first order by introducing a frame 
of reference which contains a mean orbital plzine having a constant inclination. 
Both the reference frame and the orbit plane rotate at a constant angular velocity. 
With respect to this doubly moving reference frame, the motion of the satellite 
differs from pure elliptic motions by only periodic perturbations. 

King-Hele [Ref. 26] advanced the approach taken by Roberson. He intro- 
duced a non-rotating reference frame with a orbit plane that continually rotates 
at a non-constant rate about the Earth’s axis. A relation is found between the 
rotating orbit plane and the angular rate of travel of the satellite. The equations 
of motion are written in position coordinates, but are subsequently rearranged 
in terms of a modified set of orbital elements. The inclination of the rotating 
reference plane is held strictly constant in the analysis. King-Hele formulated 
the problem by employing a power series expansion in terms of the eccentricity; 
therefore, the method is limited to small eccentricities. The final solution con- 
tains an incomplete set of workable elements and hence the method gives mostly 
a qualitative description of the satellite’s behavior due to oblateness. 

King-Hele’s analysis was the inspiration for the work of Brenner and Latta 
[Ref. 27|. They improved King-Hele’s original analysis by abandoning the condi- 
tion of constant orbit inclination and by retaining the eccentricity in closed form 
expressions. Brenner and Latta limited their analysis to small eccentricity al- 
though the method is not so restricted. They obtained an approximate first order 
solution and demonstrated that the method is valid for higher order analysis. 
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An advantage of the direct method is that one may use ordinary perturbation 
analysis (the Method of Strained Coordinates) to solve the equations of motion. In 
addition the chosen set of orbital elements xised throughout the analysis correspond 
closely with the classical elements; therefore, physical interpretation is facilitated. 

The disadvantage in using the method is shared by most all other procedures 
that must use a perturbation scheme. The process requires the manipulation of 
massive algebraic expressions. However, this drawback has been greatly reduced 
by the introduction recently of large symbolic mathematics programs, such as 
MACSYMA, that handle the bookkeeping. 
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III. THE TWO-BODY PROBLEM 



A. INTRODUCTION 

The central problem of celestial mechanics is the two-body (or Kepler) prob- 
lem. Simply stated, the problem is to solve for the motion of two particles inter- 
acting through their mutual gravitation in an isolated space. A solution of the two 
body problem often represents physical reality in an zicceptable way. For instance, 
the orbit of the Earth around the Sun may be treated, to a first approximation, as 
a two-body problem because the influence of perturbing bodies, the Moon, Jupiter, 
etc., are small compared with the Sun’s gravitational attraction. Likewise, more 
complex problems such as a spacecraft mission to Mars, a four body problem - 
Earth, Sun, Mars, spacecraft, may be treated by breaking the flight into three 
two-body problems. Such a technique is used in the patched conic method where 
two-body solutions are literally patched together. 

There are other more important reasons for studying the two-body problem. 
It is the only gravitational problem in dynamics, other than very specialized cases 
in the three-body problem, for which there is a complete and general solution, 
and it is possible to gain considerable insight into more general phenomena of 
motion by a thorough study of the two-body problem. In fact, the most complete 
theories of celestial motion use functions appearing in the solution of the two- 
body (elliptic case) problem as elementary functions. The solution is central to 
this present analysis because it will serve as a starting point for the generating of 
analytical solutions that are valid to higher orders of accuracy. These solutions, 
called general perturbation theories, are the subject of Chaper VI. 
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In this chapter, the equations for the two-body problem will be derived and 
solved. The first step is to chose a coordinate system in which the laws of Newton 
hold (an inertial coordinate system). In practice, the reference frame of the “fixed” 
stars provides a very good approximation to an inertial reference frame. Next, 
following the method of Nelson and Loft (Ref.28 pp. 82-84] it will be shown that 
the center of mass of two bodies in this coordinate system is also an inertial 
reference point. Then the equations of motion will be derived and solved using 
the polar angle 6 as the independent variable. 

B. THE DIFFERENTIAL EQUATION 

Figure 3.1 shows two mass centers at position rj and T 2 . The origin O is 
defined to be an inertial reference point. The distance between the two mass 
centers is r where 

r = F2 - Fi. 

Combining Newton’s third law of motion and his law of universal gravitation gives 

K 




Figure 3.1: Two Body Problem 
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the force equations for the two bodies: 



d^Ti 


Gmi m2 T 


dP 


r 5 


d?T2 


—Gmim2 r 


^^~dF 





where G is the universal gravitational constant. 
Adding (3.1) and (3.2) and integrating results in 



(3.1) 

(3.2) 



dvi d.T2 

mi— — hm 2 — ;— = constant. 
(O' cLt^ 



(3.3) 



Now the center of mass of the system (barycenter) in Fig. 3.1 is defined as: 



(mi + m2)ro = miFi + m2F2. (3.4) 

By combining equations (3.3) and (3.4) the following result is obtained, 

dro mi dri m2 dr^ 

= — I — = constant. 

dt mi + m 2 dt mi + m 2 dt 

So the barycenter moves with constant velocity, and it too is an inertial reference 
point. Subtracting (3.1) from (3.2) results in 



<Pt 2 d^Ti (Pr —G{mi + m^) r 



(3.5) 



dt^ dt^ dt^ f3 

the solution of which gives the position of either body relative to the other. Now 
choose the barycenter as the origin and define the position of mi and m 2 , respec- 
tively, as 

—m 2 



roi = Ti - To = 



mi -f m 2 



ro2 = F2 - To = 



mi 

] 

mi -f m 2 
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Substitution of the above equations in (3.1) and (3.2) yield septate equations for 
the motion of each body relative to the baxycenter: 



d*ro d^roi _ d^roi 

dt^ dt^ ~ dt^ 



—G m* 

(mi + 



(3.6) 



d^Tp dpTp2 _ d^rp2 

dt^ dt^ dt^ 

Equation (3.5) eind equations (3.6) are 
stant, so that 



-Gmf 

(mi +m2)^r^2 

of identical form, differing only by a con- 



dor GM r 

^ ® 



(3.7) 



is the vector differential equation of motion for either of the two bodies, r is the 
distance to the other body, or to the barycenter, according to the appropriate 
choice of GM. For the problem of a satellite orbiting the earth, the mass of the 
satellite can be neglected in comparison with the mass of the earth, therefore GM 
is the product of the universal constant and the mass of the earth. 



C. THE INTEGRATION OF THE TWO-BODY PROBLEM 

There are no cross products involved in equation (3.7); therefore, all motion 
must lie in the plane that contains r and The scalar components of acceleration 
are 



r 



dt^ 



+ 



2drde 

dt dt 



(3.8) 



dV Uey_-GM 



(3.9) 
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Writing (3,8) as = 0 and integrating yields 

<tdO , , , 

r — = h = constant (3.10) 

at 

where h is the specific angular momentum. 

Equation (3.10) is an exact integral of (3.8). It corresponds to Kepler’s 
empirical law of constant areal velocity which states that the area swept out by 
the radius vector of a planet is uniform in time. 

From (3.10), the independent variable may be chzinged from t to 6, e.g., 

Ct hr dr 

dt dO 

so (3.9) becomes 

1 d^r 2 /dry 1 _ GM __d? 1 _ GM 

Since this equation is linear in the reciprocal of r, it may be written as 

d}u GM 

which has the solution 

u — + i4cos(^ — w) (3-11) 

or reintroducing r, (3.11) becomes 

h^GM 

1 + [Ah'^ /GM) cos{0 — ly) 

which may be written as the polar equation of a conic section 
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(3.12) 



so that 



P 

1 + e cos(5 — w) 



p = h^GM and c = Ah^GM. 

D. ELLIPTIC MOTION 

Equation (3.12) is the equation of a conic with prime focus at O. The conic 
has a semi-latus rectum h^/GM, an eccentricity e, and a semi- major axis a that 
makes an angle f = 6 — u with the horizontal axis (Figure 3.2). The extreme 
endpoints of the major axis of the orbit are referred to as apsides or apses. The 
point nearest the prime focus is called perigee aoid is given by 6 = u>. The point 
farthest from the prime focus is given by ^ = w -t- 180® and is called apogee. The 
angle w, “argument of perigee”, will be discussed later in this chapter. 

The energy of the satellite in the orbit is conserved and is equal to 

E = m,(v^/2 — GMjr) = m,C 




Figure 3.2: The Elliptic Orbit 
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where C is the total energy of the satellite, v^/2 the kinetic energy ajid —GMfr 
the potential energy of the satellite, all per unit mass (m,). In general C is equal 
to —GMf2a, so that the satellite’s energy depends only on the semi-major axis. 
From this relationship it is easily shown that the satellite’s velocity is given by 

= GM(2/r- 1/a). 

Then the relationships Vp = a(l — e) and = a(l + e) result in 

, GM(1 + e) 

^ a(l - c) 

and 

, . GM(I - e) 
a(l+e) 

so that the velocity is a maximum at perigee and a minimum at apogee. 

The area of the ellipse is — e^), and the rate of description of area is 

= |. Since = GMa{l — e^), the orbital period may be written as 

T = 

By defining the mean motion n as T = ^, so that n^a® = GM, one may proceed 
to derive an expression for position versus time in the elliptic orbit. 

The orbital ellipse APB with center at C touches at perigee. A, and at 
apogee, B, a concentric circle also centered at C which has as a radius the semi- 
major axis of the ellipse. The circle C is known as the auxiliary circle, and is 
geometrically related to the ellipse by the relation 

PN = P'N^{1 - e2) 
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where c is the eccentricity, P any point on the ellipse, N the foot of the perpendic- 
ular through P upon AB, and P' the intersection of this perpendicular with the 
circle C. The angle ACP' = ^ is known as the eccentric anomaly (Figure 3.3). 

Let T be the time of perigee passage and t the time, then t — t is the time 
since perigee passage. The quantity. 



M = n(t-r), (3.13) 

is called the mean anomaly. Using the geometry of Figure 3.3, one may derive the 
following relationship between the anomalies: 

M = E — esinE, (3-14) 

which is known as Kepler’s equation. If the position of a satellite in a fixed orbit 
relative to the earth is desired at some specified time t, then equation (3.13) gives 
M and (3.14) gives E . The distance to the satellite is found by the relationship 

r = a(l — ecos E). 

The angular position of the satellite is defined by the true anomaly, 6 — w — /, or 
the angular distance from perigee: 




Figure 3.3: The Eccentric Anomaly 
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tan(//2) = (^) tan|. 



In actual practice, the second step in the above process, solving (3.14) is a bit 
more involved since a closed form solution to Kepler’s equation does not exist. 
However, dozens of methods of successive approximations have been devised. For 
instance, by use of a numerical method M can be calculated for a few values of E 
and then the correct E can be found by inverse interpolation. 

E. CONSTANTS OF THE MOTION 

In section B, the original equations of motion (3.1) and (3.2) were reduced 
to (3.7). Thus the problem was reduced from one of three second order differential 
equations requiring twelve constants to one of three second order equations with six 
constants. A discussion of the constants of motion, some of which were introduced 
during the solution of the two-body equation, is the purpose of this section. The 
six constants may be written in a variety of forms, a choice among forms is usually 
made with the purpose of simplifying the problem. 

Equation (3.7) was solved by the classic technique of changing the indepen- 
dent variable from t to 6. But (3.7) in its present form is integrable; that is, 
there exists sufficient time independent first integrals, or functions that are con- 
stant along the motion, to specify each orbit. The first of these is the angular 
momentum. Cross-multiplying (3.7) by r results in 



or 



dV GM 

■•X *!+■•>= —■■ = 0 



d^r 

r X = 0. 
dt^ 
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And since 



d ( dr\ dr dr <Pr 

— rx— = — X — + rx -r-r 

dt \ dtj dt dt dt^ 

then 

|(r^v) = 0. 

By integration 

r X V = h 

where h, the angular momentum, provides three constants of motion. Similarly, 
by cross-multiplying (3.7) by h, one may obtain the Lenz vector 

dr h T f de \ 

® ^ GM ~r “ j 

where e is a vector along the major axis of the orbit pointing toward the position of 
closest approach or perigee (|e| = e). The vector e provides only two independent 
constants since h and e are perpendicular vectors (h • e = 0), and one remaining 
constant is required. 

The vector integrals h and e specify that the orbit will lie in the plane per- 
pendicular to h and have a shape determined by e. The classical orbital elements: 
a (semi-major axis), e (eccentricity), i (inclination), fl (longitude of the ascending 
node), and u (argument of perigee), may be derived from these vectors. 

From (3.12) the equation for r may be written as 

_ h^GM 
1 -f e cos / 

where / is the angle between e and r (the true anomaly). 

Restricting c to the elliptic case results in 
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r\ = 






<r < 



GM[l + e) 
so that the major axis of the ellipse is 

2a = T\ "i" T2 — 



GM{1 - c) 



= ^2 



2h} 

GM(1 - c2) ■ 



Knowing e and a gives the shape and size of the orbit. The orbit now can be 
oriented in a coordinate system. Reference is made to Figure 3.4 and the two 
angles i and fl. i is the inclination of the orbit plane defined as the angle between 
the equatorial plane and the orbit plane. Since this is the same angle as the angle 
between the z axis (k unit vector) and the angular momentum vector h, i may be 
found by 

. h, 
cos I = — . 
h 



k 




\ 



Figure 3.4: Constants of the Orbit 
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The angle Q, the longitude of the ascending node, b the angle in the equato- 
rial plane, between the i unit vector and the point N where the satellite crosses 
through the equatorial plane in a northerly direction measured counterclockwise 
when viewed from the north side of the equatorial plane. The vector n lies along 
NO such that 

n = k X h. 



Therefore fl may be found by 

r> 

cos n = — , 
n 

As was stated above, the vector e, points toward perigee. The angle u (the 
argument of perigee) measures the distance between NO and perigee and can be 
found by 

n • e 

cos u = . 

ne 

With the constants a,e,i,Cl,u specified, the orbit is defined in the coordinate 
system. The remaining task, to locate the satellite in the orbit at any time, 
requires one more constant. 

The final constant of motion is given by the relationship between the mag- 
nitude of the angular momentum and the true anomaly (/): 

/ hdt = f r^df. 

Jt •'0 

By a change of variables from / to the eccentric anomedy E, the above equation 
may be easily integrated. The result by this analytical method is identical to 
equation (3.14), Kepler’s equation, which was previously derived by geometry. 
The constant t, time of perigee passage, is the final constant of integration. 
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F. SUMMARY 



This completes the analysis of the two-body problem. It has been shown that 
a combination of six constants will strictly define the motion of a satellite under 
the influence of a central gravitational force. The six that were chosen, referred 
to now as the orbital elements, can be \ised to find other constants of the motion 
including r and v, the distance and velocity vectors of the satellite. 
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IV. SATELLITE PERTURBATIONS 



A. mTRODUCTION 

As was demonstrated in the last chapter, the classical two-body problem has 
solutions that can be written in closed form when the polar angle (or the eccen- 
tric or true anomaly) is used as the independent variable. If an additional force 
acting on either of the two bodies is introduced, the resulting equations of motion 
usually no longer have closed-form solutions. When the magnitude of such a force 
is small compared to the central gravity term, the force is called a perturbation. 
The resulting orbit experiences a small departure from the Keplerian orbit, at 
least initially. These departures are also called perturbations. Under certain cir- 
cumstances, it is possible to make analytic approximations to the effects of the 
perturbing forces, though a precise solution cannot be obtained. Generally, the 
methods consist in determining the exact equations of motion and then assuming 
that their solutions do not depart appreciably from the case of no disturbing force. 
Then only an indication of the actual motion of the body can be obtained. Pre- 
cise solutions can be found for specific initial conditions by numerical integration 
techniques, but these solutions give little insight into the dependence of the mo- 
tion on the parameters of the disturbing force. In some cases, the approximations 
obtained with analytic methods may exceed the precision of numerical methods, 
especially if the prediction is required for a long period of time and there is a clear 
dominance of one particular force. 
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In the case of a close satellite about a non-spherical planet, a potential func- 
tion V can be formed such that 



V = Vo + R 

where Vq is the potential function due to the two-body problem and R the dis- 
turbing ftmction that is at least an order of magnitude less than Vq. Many general 
perturbation theories make iise of the fact that the two-body orbit of the body due 
to Vo changes slowly due to R, auid they attempt to obtain analytical expressions 
for the changes in the orbital elements due to R within a specified time interval. 
If the elements of an elliptical orbit are OojCo, and ^o> the ellipse with 

these elements is referred to as the osculating elements at to. The velocity of the 
disturbed satellite at this time in its osculating ellipse is equal to its velocity in 
the actual orbit. 

At a future time ti, the elements will change due to the presence of R. 
For instance, Oq will be changed to uj, and the quantities (uj — ao)»(ci — Cq), 
etc. are called perturbations in the elements. Corresponding to the perturbations 
in the elements are perturbations in the coordinates and velocity components. 
There are, however, at least two reasons for using orbital elements rather than 
coordinates to describe the motion of the satellite. First, the elements do not 
exhibit a variability of anomalistic motion that the coordinates do, therefore any 
variation can be attributed directly to the perturbing forces. Second, the elements 
give a clearer geometric picture than do the coordinates; hence the effect of the 
perturbation on the orbit can be seen immediately. 

There are various kinds of disturbances that an orbit can experience, the 
severity of each is usually due to the altitude of the satellite. It is the purpose of 
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this chapter to give a qualitative description of the most important disturbances, 
and to relate the relative magnitudes of eaxh. 

B. THE EARTH’S GRAVITATIONAL FIELD 

The two-body problem assumes that the earth is a sphere, while in reality the 
earth is flattened somewhat at the poles and bulges correspondingly at the equator. 
Such a shape is called an oblate spheroid. In the science of geodesy, it has been 
useful to define a reference ellipsoid as a mathematical surface which is an idealized 
approximation to the earth’s actual surface. The study of satellite orbits has 
established a flattening of the terrestrial ellipsoid as 1/298.24, which corresponds 
to a difference between the equatorial and polar radii of 21.4 kilometers. 

Another surface commonly used in geodesy is the geoid, which is the equipo- 
tential surface that coincides on the average with mean sea level in the oceans. The 
geoid is everywhere perpendicular to a plumb-line since gravity is always normal to 
its surface. Before the advent of artificial satellites, it was generally accepted that 
except for relatively small variations resulting from the presence of mountains or 
deep valleys, the geoid could be regarded as approximately an ellipsoid. Now, the 
surface of constant gravitation can be more accurately portrayed by representing 
the potential as a series of quantities known as spherical harmonics, each of which 
makes a contribution, positive, negative, or zero, to the total. The contribution 
of any harmonic is determined by a factor, usually represented by the symbol J 
and commonly referred to as the value of that harmonic. These J values for a 
planet’s gravitational field can be determined from observations of a satellite orbit 
and they can be related to the shape of the geoid. 

A large number of harmonics may be required to precisely represent a planet’s 
gravitational field, but in practice the higher harmonics make such a small contribution 
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that they can be neglected, at least to a first approximation. The zonal harmonic 
Jo expresses the overall size of the geoid, while Jj, the first degree harmonic de- 
termines the center point of the geoid in the north-south direction. The other 
harmonics represent deviations from the spherical shape as shown by Figure 4.1. 
It is seen that the contributions from the even harmonics are symmetrical about 
the equator, while the odd harmonics corresponds to contributions that are asym- 
metrical. The degree of the harmonic gives the number of undulations in the shape 
of the surface. 




Figure 4.1: Qualitative Representation Of The Harmonics Of The Geoid 

The results given thus far consider only the zonal harmonics, which are in- 
dependent of longitude. The tesseral harmonics give east-west deviations from 
symmetry. Satellite observations of the tesseral harmonics have led to the con- 
clusion that the equator of the earth’s geoid is slightly elliptical, the difference 
between the longest and shortest ajces being about 400 meters. 

The following tables from Kozai [Ref. 29] gives representative values for the 
coefficients of the earth’s zonal and tesseral harmonics. 
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TABLE 4.1: ZONAL HARMONICS 



n Jn X 10^ n Jn X 10^ 

2 1082.48 ± i04 6 .39 ± i09 

3 -2.57 ± .01 7 -.47 ± .02 

4 -1.84 ± .09 8 -.02 ± .07 

5 -.06 ± .02 9 .11 ± .03 



TABLE 4.2: TESSERAL HARMONICS 




C. ATMOSPHERIC DRAG 

Atmospheric drag is the most complex and the most difficult of the artificial 
satellite perturbations. The complexity arises from the fact that an exact force 
law is not known and the atmosphere is variable. This variability results from 
the fact that the atmosphere is not spherically symmetric and that lunisolar tides, 
diurnal heating, strength of the solar wind, etc., all effect atmospheric density. 
In addition, the atmosphere is moving with the rotating earth. Therefore this 
perturbation is difficult to model analytically or numerically. 

For most Earth satellites, atmospheric drag removes the satellite’s energy, 
thereby causing it to drop in altitude and thus increase its speed. The increased 
speed and the lower altitude increases the drag force with the result that the 
satellite spirals into Earth. 
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Satellite drag is measured along the same direction as the velocity vector of 
the satellite. The formula used is 

Fj 3 = - v.|(v. - v) 

where v is the satellite’s inertial velocity, v* is the inertial velocity of the atmo- 
sphere, p is the atmospheric density, A is the effective cross-sectional area of the 
satellite, and Cd is the drag coefficient, Cp « 2. 

D. SOLAR RADIATION PRESSURE 

In contrast to drag which is a tangential force, solar radiation pressure is 
radial. By solar radiation on a satellite is meant the net acceleration resulting 
from the interaction (i.e., absorption, reflection, or diffusion) of the sunlight with 
the surface of the satellite. In general, the illumination of a low altitude satellite 
orbit will be non-uniform due to the Earth’s shadow, albedo, and atmospheric 
absorption. The magnitude of the solar radiation pressure is given by 

AP 

4nmcD‘^ 

where A is the effective cross-sectional area of the satellite, P is the total radiated 
solar power, m is the satellite’s mass, c is the speed of light, and D is the satellite- 
Sun distance. 

It can be seen that the perturbing effect of solar radiation particularly effects 
satellites with a large area to mass ratio. Such was the case of the satellite Echo I, 
which when inflated was a sphere with a total exposed area of about 31,000 square 
feet and weighed only 135 pounds. Its initial orbit was approximately circular with 
an altitude of 1000 miles, but the pressure of solar radiation brought down the 
perigee to 600 miles at times. 
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A. OTHER PERTURBING FORCES 



Gravitational perttirbations due to other celestial bodies (the Moon, the Sun, 
and the planets) are caused by the differences between the force on the Earth and 
that on the satellite or the “tidal” force. Perturbing forces from the Sun and Moon 
become significant as the altitude of the satellite increases. Planetary perturba- 
tions are very small, with Venus and Jupiter providing the largest contributions. 

Since classical mechanics is used in this analysis, relativistic corrections may 
be regarded as a small perturbation to the motion. Other minor perturbing forces 
include atmospheric lift and electromagnetic forces. 

The following table from Milan! [Ref 30] gives the magnitudes of various 
perturbations on three different satellites; 

TABLE 4.3: ACCELERATIONS ON SPACECRAFT AT VARIOUS 

ALTITUDES (cm/sec*) 



Cause 


G eosynchronous 
satellite 

a = 42,160 KM 
A/M = .1 cm^/g 


LAGEOS 

12,270 

.007 


SEASAT 

7,100 

.2 


Earth’s 

monopole 


2.2 X 10^ 


2.8 X 10^ 


7.9 X 10^ 


Earth’s 

oblateness 


7.4 X 10-“* 


1.0 X 10"^ 


9.3 X 10-1 


Higher 

harmonics 


4.5 X 10-^° 


8.8 X 10"® 


7.0 X 10-1 


Moon 


7.3 X 10"“* 


2.1 X 10-* 


1.3 X 10-* 


Sun 


3.3 X 10““* 


9.6 X 10-® 


5.6 X 10-® 


Venus 


4.3 X 10"® 


1.3 X 10“® 


7.3 X lO-o 


General 

Relativity 


2.3 X 10-® 


9.3 X 10“® 


4.9 X 10-^ 


Air Drag 


0(?) 


3.0 X 10-10 


2.0 X 10-® 


Solar 

Radiation 


4.6 X 10"® 


3.2 X 10"^ 


9.2 X 10-® 


Earth’s 

albedo 


4.2 X 10"® 


3.4 X 10“® 


3.0 X 10-® 
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V. DEVELOPMENT OF THE EQUATIONS OF 

MOTION 



A. PRELIMINARIES 

1. Overview 

In this chapter, the groundwork will be laid for solving the equations of 
motion for a satellite about an oblate planet. To begin, a discussion of the assumi>- 
tions made in the analysis will be given. Second, the special coordinate system, 
which was introduced by King-Hele [Ref. 26] and refined by Brenner and Latta 
[Ref. 27] will be developed, followed by a derivation of the relationships among 
the astronomical angles of the coordinate system and the angles of a spherical 
coordinate system. Next, an expansion for the potential of a planet modeled as 
an oblate spheroid will be derived. Finally, the equations will be transformed so 
that they can be solved by an ordinary perturbation method. 

2. Assumptions and Limitations 

The basis for the assumptions made in this analysis is the requirement 
that a solution to the equations of motion for a satellite orbiting a planet be 
accurate to within a relative error < 10“®. Therefore, all perturbative forces 
that are of magnitude 10“® or smaller compared to 1 may be neglected. This 
requirement then allows one to model the earth as an oblate spheroid with an 
axially symmetric gravitational potential. This is to say that all zonal harmonics 
except the second in the expanded potential may be neglected. In addition, all 
coefficients of the tesseral harmonics are small enough to neglect. 
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The assumption of an axially symmetric potential is a reasonable as- 
sumption if the earth is used as an example of an oblate planet. As was suggested 
in Chapter IV, if one normalizes the expanded geopotential of the earth such that 
the contribution from the central force (l/r*) is 1, then the contribution from the 
zonal harmonic for the oblateness contribution is « 10“*, and all other harmonics 
(zonal and tesseral) are < 10“®. 

The gravitational effects of the Sun and Moon are also neglected. This 
assumption will remain valid out to a distance of at least 6,000 kilometers above 
the Earth, where the oblateness perturbation is still about 1000 times larger than 
the attraction of the Sun. However, at altitudes near geosynchronous (35,800 km) 
the perturbative forces of the Sun and Moon are nearly equal to that of oblateness, 
and therefore would have to be considered. 

The most important limitation to the analysis is the neglect of atmo- 
spheric drag. For high altitude satellites with small eccentricity, the neglect of air 
drag is unimportant; however, for very low altitude satellites the effect of air drag 
would begin to dominate the oblateness force after a number of revolutions. Also, 
high altitude satellites with highly eccentric orbits would be greatly affected by 
drag as they pass through perigee. For these particular cases, the desired accuracy 
for 1000 revolutions would not be achievable; however, in many cases the solution 
could be accurate for a few orbits. As was discussed in Chapter IV, air drag re- 
mains the most difficult perturbation to model. Since an accurate geopotential 
model allows for a better determination of fluctuations in drag of the atmosphere 
through which the satellite is moving, the analysis conducted here is valuable even 
when air drag becomes the dominating factor. 

Also neglected in this analysis, are the remaining perturbative effects 
mentioned in Chapter IV. Solar radiation pressure, the Earth’s albedo radiation 
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press\ire, and the general relativistic correction to Newton’s law of gravity are all 
neglected since their relative contributions are small. 

The final limitation to the analysis is the assumption that the differential 
equations which describe the motion of the satellite have a solution asymptotic to 
the form specified in the analysis. As long as the true solution does not depart 
from the form specified, the analysis is valid. 

3. Order of Magnitudes 

Throughout this analysis, reference will be made to the relative mag- 
nitudes of the perturbing forces; hence it is important to establish a convention 
for orders of magnitudes. The approach established by King-Hele [Ref 26] will be 
followed here with the exception that a small eccentricity is not assumed. Here, 
J « J 2 , the coefficient of the second zonal harmonic in the potential is used as the 
basic small unit. Mathematically, 

f{JJ) = 0{J^6^) 

means that there exists for all 

0 < J < 10"® and 0 < ^ < (27r)l0® 
a K independent of J and 6 such that 

1/1 < 

The following terms are defined: 

Zero order = 0(1) ~ 1.0 

First order = 0(J), O(J^0), O(J®0^), . . . ~ 10"® 

Second order = 0( J^), 0(J®^), 0( ~ 10~® 
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B. THE COORDINATE SYSTEM 



Figure 5.1 is an inertial reference system of spherical coordinates with the 
prime direction pointing toward the vernal equinox at epoch 1900.0. The equato- 
rial plane, latitude (<5), and longitude ((^), as shown in the figure have their usual 
meanings. 




Figure 5.1: The Reference System 



As was demonstrated in Chapter III, the path of a satellite of a strictly spherical 
planet lies entirely in a fixed plane, and the motion of the satellite is described by 
the solution to the two-body problem. With the angular momentum vector fixed 
in space, the intersection of the orbit plane and the equatorial plane describe the 
fixed angle fl measured from the prime direction to the intersection (or node). 
The longitude of the ascending node (fl), the inclination (i), and the argument of 
perigee (w), fix the orbit plane in the coordinate system. 
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The effect of the oblateness perturbation is, in general, to move the satellite 
out of the original two-body orbital plane. One specific effect is to rotate the 
orbit plane about the planet in the opposite direction to the satellite’s motion. 
The physics of this motion is easily demonstrated in Figure 5.2. Oblateness is 
represented by an additional mass about the equator of the planet. The radius 
vector r and the satellite S lie in a plane which is named here the “reference 
plane”. If there were no equatorial bulge, the direction of the gravitational force 
would coincide with the radius vector and the angular momentum vector would 
remain in a constant direction normal to the plane. The equatorial bulge; however, 
adds a component of force that does not lie along r. This additional force adds a 
torque r = r x F. The direction of r is into the paper at S for a prograde orbit. 
Therefore as the satellite orbits, the angular momentum vector rotates about the 
Z <ixis. The reference plane also rotates, and the rotation is measured by the 
change in fl. 

Polor 

axis 




Figure 5.2: Rotation of the Reference Plane 



The line NON' in Figure 5.1 is referred to as the “line of nodes”, and the 
ascending node describes the point where the satellite enters the northern hemi- 
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sphere. For the two-body problem with a non-rotating plane, the line of nodes is 
a fixed reference line. Now, for the rotating plane, when the satellite is actually 
at a node, the line of nodes passes through the satellite, but at other times the 
definition “line of nodes” is an arbitrary one, and the angle fl is simply given as 
a function of time. 

Another motion of the reference plane that is a result of the oblateness 
perturbation is an in-plane rotation called precession. As was demonstrated in 
Chapter III, a particle executing bounded, non-circular motion in a central force 
field will always have a radial distance from the force center that is bounded by 
’’mor ^ ^ »’mm- That is, r is bounded by the apsidal distance. Such an orbit is 

called a closed orbit and it is characterized by the fact that all apsidal angles are 
equal. For instance, the apsidal angle is tt for elliptic motion. But if the radial 
dependence deviates slightly from l/r*, then the apsides will precess or rotate 
slowly in the plane of motion (Figure 5.3). This motion is analogous to the slow 
rotation of elliptic motion of a two dimensional harmonic oscillator whose natural 
frequencies for each dimension are almost equal. The rotation of the line of apsides 
is also known as the precession of perigee since the value for 6 for which r is a 
minimum varies in time as the apsides rotate. 




Figure 5.3: Rotation of the Line of Apsides 
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For a satellite orbiting an oblate planet the rotation of the line of apsides is 
opposite the direction of satellite motion for inclinations less than sin“' ^4/5 and 
in the same direction as the satellite for inclinations greater than sin' If 

the inclination is exactly sin“^ y^4/5 then there is no rotation of the line of apsides, 
and perigee remains constant. This inclination is known in Astrodynamics as the 
“critical inclination” . A detailed discussion of the critical inclination is contained 
in Appendix C. 

A final motion of the reference plane caused by the oblateness perturbation 
is the periodic variation of the inclination about the initial inclination t’o- 

C. ANGLE RELATIONS 

As a preliminary to deriving the equations of motion, it is necessary that 
the relationships among the spherical coordinates and the angles describing the 
reference plane be established. 

Reference is made to Figure 5.4, where and j' define the equatorial 

plane. Let a, b, and c denote three unit vectors, where b and c are in the reference 
plane, c points to the initial point 0 = tt/2 in the reference plane where 0 is 
measured from the line of nodes, b points to the position of the satellite, S. a is 
perpendicular to the reference plane and therefore points in the same direction as 
the angular momentum vector h. a and c are both perpendicular to the line of 
nodes. Then 



b = cos <f> cos 6 i + sin^ cos^ j + sin 6 k. 

Now measure b from the line of nodes. 

b = cos(<j!> — n) cos £i' -f sin(^ — U) cos^j' + sinik. 
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a and c are therefore 



a = — sin tj* + cos *k 

c = cos » j' + sin t k. 

Portion of the Reference Plane 




Since a and b are perpendicular, a • b = 0, therefore 

tan^ = tant sin(^ — fl). 



Continuing 



results in 



b • c = |b| |c| cos(^ — 7 t/ 2) 
sin((^ — n) cos ^ cosx + sin^ sini = sin0. 



(5.1) 
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And using the relation from equation (5.1) results in 



sin^ = sin^ sint. (5.2) 

Finally, from c x b, the following relationship is obtained 

cos 6 = cos 6 sec (» - n). (5.3) 

(5.1), (5.2), and (5.3) are the required angle relations. 

D. THE POTENTIAL 

In Chapter III, a simplified gravitational potential GMjr weis used with the 
assumption that the attracting body was spherically symmetric. This simplified 
potential caused the satellite to move in conic orbits. As has been stated, the 
planets are not spherically symmetric but are bulged at the equator, flattened at 
the poles, and are generally asymmetric. An expanded expression for the gravi- 
tational potential will be developed in this section. The final expression for the 
potential is subject to the assumption made in Section B of this chapter that the 
planet about which the satellite revolves is approximately an oblate spheroid. 

Now, regardless of the nature of an attracting body, the potential V must 
satisfy one of the following differential equations. For regions within the attracting 
matter 

= Ai^pG where p is the density (5-4) 

which is Poisson’s equation. 

For regions outside attracting matter 

VV = 0 (5.5) 



which is Laplace’s equation. 
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This discussion will be restricted to Laplace’s equation which, if V is written 



as a function of the spherical coordinates may be written as 



dr y dr / **” r* cos 5 d5 \ d5 y **” r* cos* 6 d<f>^ 

Any conservative force field F can be written as the gradient of a potential. 
Equations (5.5) and (5.6) state that div (grad V) is zero. The assumption that F 
is rotationally symmetric mezuis that dV/d<f> is zero. Therefore a solution of (5.6) 
is sought that is a product of a function of r alone and a function of 6 alone: 




V{r, 6 )=h{r)-M 6 ). 



By multiplying (5.6) by r*, dividing by fi{r)f 2 { 6 ), and rearranging terms, the 
following result is obtained. 



1 d 

h dr ' ’■* 



dr 






/j COS 6 dS 

Since the left side is a function of r alone and the right side is a function of 6 
alone, both members are equal to a constant, say n(n + 1). fi must satisfy 



which has the solution 



/n = Ar-"-^ + Br". 



It is desired that this function be zero at r = oo, and that it be analytic and 
single-valued there. Therefore n is chosen to be an integer greater than zero, and 
B most equal zero. Therefore 

fi=Ar—\ 



The equation for /j is 



A. 

dS 




-h n(n + l) cos 6 / 2=0 
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which may be rewritten as 



dsinS 



(1 — sin^5)— = 



df2 



d sin S 



+ n(n + l )/2 = 0. 



(5.7) 



This is known as Legendre’s equation. Using the method of Frobenius, a series 
solution to Legrendre’s equation is found to be 

/2 = CP„(sinS). 

Therefore the solution to (5.6) may finally be written as 

V(r,^) = (sin 5) 



where C is chosen as 1. 

There are many ways to consider the polynomials Rodrique’s formula for 
P„(sin(5) is 

P„(sin5) = 7 ;;^ j,.!Tcx„ (sin^^-l)’’. 



so that 



2"n! d(sin^)" 



Po = 1) Pi=sin(5, P 2 = -(3sin^ (5 — l), etc. 



The general solution of (5.6) is then written as 

^(^>^) = X^>inr""~^P„(sin<5). 

n=0 



(5.8) 



The above mathematical derivation of an expression for the potential may 
be enhanced by a more physical formulation. Referring to Figure 5.5, let a unit 
mass m be placed at a point P which is a distance r from the center of mass of 
a bounded distribution of total mass M. Let dm be an element of the mass at a 
distance t from O. Then the potential at P due to dm is 

-K^dm 



dV = 



(£2 + — 2r£ cos 6) 
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and the total potential is 




Figure 5.5: Potential at an Exterior Point Due to an Irregular Mass 

The denominator of the right side of (5.9) may be expanded such that 

(£^ + r^-2r£cos^)-^/^ = if; Pn{sin6) 

where 6 = tt/2 — 6. 

Therefore the following may be written 

V = — ^ [ dm — f £ sin 5dm + [ (}dm 

r Jm Jm 2r^ Jm 

f s2 . 2c j 

— I Sin 6 dm + 

2 Jm 

A description of the physical significance of each of the above integrals can be 
made. The first integral is the total mass, the second is the first moment about 
an axis through O perpendicular to OP and is zero when the origin is chosen as 
the center of mass. The third integral is the moment of inertia about the origin, 
and the last integral may be written as 




where Iq is the moment of inertia about the origin and I is the moment of inertia 
about the line OP. 
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The potential is then 



K^M _ ^ 

r 2r» 



(2/0-3/). 



Now the first term is the potential due to a homogeneous solid sphere. The second 
term arises from the departure of the mass M from spherical shape. 



With the physical description established, (5.8) may be written as 




(5.10) 



Note that only Pq and Pj remain. 

An advantage of (5.10) is that it can immediately be written down once axial 



satellite about the Earth) can be used to determine (or ^3,^4, etc.; for general 
potential). 

E. THE EQUATIONS OF MOTION 

Referring again to Figure 5.1, the components of the velocity of a satellite in 
the coordinate system are: 



symmetry has been assumed, and any suitable experiment (such as an orbit of a 



V = Vr + + Vs 



or 




An expression for the kinetic energy {T) is 




( 5 . 11 ) 



The equations of motion may be written using Lagrange’s equation 
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dd[T-V) alT-V) „ , ^ 

-—^-—^ = 0 q = T,6,oTd> (5.12) 

where q = ^,T \s the kinetic energy eind V is the potential energy. 

Applying (5.12) to the expression for T (5.11) and V (5.10) results in 



dV {dSy 2c(d<i>Y 

— r — I — r cos <5 — = - 

dt^ \dt J \dt J 

d ( ^ds\ j . , , (d<t>y 

^fr=cos"5^) =0 

dt\ dt) 



dV 

dr 

dV 

dS 



(5.13) 

(5.14) 

(5.15) 



The last of these equations gives an integral which can be used to eliminate 
time from the system, i.e., 



2 2 r d<t> y 

r COS b—— = h cos »o. 
dt 



Let r = p/u then 



d h cos t’ou* sec^S d 

Jt ^ 



p2 d4> 

Starting with equation (5.14) and applying (5.16) to the first term results in 



(5.16) 



d ^ 2 dS\ _ d ( r^h cos^ j‘oscc*(5\ d6 

dt \ dt j dt ^ J dt^’ 

and continuing to apply (5.16) so that (5.14) finally becomes 



or 



d^ tan S dV r* cos* S 

d(f>^ dS cos^ t'o 



d* tan 5 
d4>^ 



+ tan 6 — 



sin 8 cos^ 8 
cos^ I'o 



(2 Ju) 



(5.17) 
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where 




p = h^/GM. 

It is desired to get the left side of (5.17) in terms of the independent variable 
6. To begin, the right side of the angle relationship (5.1), 



tan 6 = tan i sin(<^ — fl). 



is substituted into (5.17). 
Then using (5.3) 



cos 6 = cos 6 sec{4> — fi) , 



(5.17) becomes 




where primes denote differentiation with respect to 6. 

Later, <f> will be eliminated from (5.18), but first the remaining equation of 
motion (5.13) will be rewritten in a form like (5.18). 

To rewrite equation (5.13) in a useful form, the independent variable is again 
changed by use of (5.16). Multiplying the expression by and noting that 

(tan^ ^)' = 2 tan S sec^ 

du 
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results in: 



u 



M 



+ 



u 



<j)'^ cos* 6 + 



(tan 6)'^ 
(sec^ 6) 



^ u'(tan* Sy 

sec^ 6 



u'<f>" 

<!>' 



<f>'cos*6dV 
h}u} dr ■ 

Evaluating ^ and multiplying (5.19) by results in: 



(5.19) 



u" + u 



cos* t 



COS'' to 



— [l + J u*(l — 3 sin* 6)] + u" 



(l + ~ *' ^ ^ *) 



u' cos* i sec^ S 






2 sin 6 cos 6 (tan 6)' 



4 >'. 



— u 



cos* t sec* 6 + 



(tan ^)'* cos* i 



- 1 



(5.20) 



Equations (5.18) and (5.20) are the equations of motion. Before proceeding 
to solve these equations, it is necessary to remove <j>. Combining equations (5.1) - 
(5.3) results in the following expression: 



tan(^ — n) = cos t tan 6. 



Differentiating this expression with respect to 6 results in 



4 >' = 



cos t 
cos^ S 



+ n'- 



t' sin 6 cos 6 sin t 

cos^ 6 



(5.21) 



or 



J_ _ cos*^ 1 

4>' cos t 1 + Q<£ 02 il _ sin 6 cos B tan i 

COS 1 

From (5.2), sin*^ = (sin 6 sin t)* or 

cos*^ = 1 — (sin B sin t)*. 



(5.22) 
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Substitution of this expression plus (5.21) and (5.22) in (5.18) and (5.20) 
results in completely general equations of motion in terms of the dependent vari- 
ables: 1, fl,y; and the independent variable 6. 



® * (^) ^ 3 sin t ^ ^ cos 0 sin^ 6 



cos 1 



cos t 



+ 

-I- 



3 sinS'^g^cos^^ sin 0 6 sin’t cos^ 0 sin 0 

COS t cos X 

3 sin* ^ cos* 0 sin 0 2 sin* ^ sin 0 

cos t cos i 



■ di d*n . , 

-h 3 cos 1— -77:7 sin 0 
d0 d0^ 

. 2 



ry di _• n 
^d0 d$^ ° 

cos t 



3 sin*t (^)%in 0 . _ . f dClV . ^ 

-I- — 7 1- 6 cos t sin 1 I — 1 sin 0 



cos t 



3 sin t (S)%in0 ^ sin* i§k§5m0 



COS X 



d$ 

COS t 



.d*t dn . , S^sin0 
+ 3cos t— — sin 



d0^ d0 



— 2 sin t— sin 0 
cos t d0 



d^i . 6 sin sin 0 3 sin* cos* 

d0'^ ^ cos i cos t 



12 sin* cos* 0 . . .dn d*n 

_l_ 3 cos t sin t— — — 

cos i d0 d0^ 

. ,d*n ^ 2sin*t^^cos0 

+ sint-— COS0 + — 



cos 0 



d0^ 

. di dfl 

2 cos* i sin iJu sin 0 
cos^ I’o 



cos t 
cos t 



— 2-— cos 0 
d0 



(5.23) 



and equation (5.20) is: 
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fu _ 3 sin^‘ (g)^^sin^(2g) 

d6"^ 4 cos^ i 

3 sin(2^) 

cos^t 

sin(2^) 

2 cos ’ i 

2 sin* cos 0 sin 0 sin (2^) 
cos’ i 

sin* i u cos 6 sin B sin(2z) sin* i ^ cos' 
cos ’ i 2 cos’ t 



<^os*^sin(2^) 3 sin sin(2^) 



+ 



cos’ i 



sin tSgTsin{2^) , sin sin(2^) 



+ 



d\ d*n du , 



de ds ds^ 
cos’ i 



• de dS'^ 
cos i 



+ 



2 cos’ t 

3 sin^ t ( 4i“ sin^ 6 ■? c;r,< , <*0 <l’n du _• 4 a 

stn t sin u 3 sm t^^^sin 6 



cos’ i 



cos’ X 



cos’ I cos ’ i 

2 Sin CQs sin ^ 

cos’ t 

6sin'*:(^^ gcos^sin*^ 
cos’ i 

^ sin*t0g^cos e sin*g 
cos’ t 

^sin^^(g)^ geos g sin* g 
cos’ t 

9 sin* : g g g cos’ 6 sin* ^ 
cos’ i 

2 sin* tggjg cos’ ^ sin*^ sin*ig^u cos’^ sin*^ 
cos’ t cos’ : 



^sin(2z) 
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+ 

+ 

+ 



6 sin^ * (d») 2 sin^t#^sin^0 



cos* I 



dB dB‘» 
COS i 






' de de^ de 
cos^ i 



cos t 



2sin^*gf^sin^0 . 2 sin f sin^fl 



* de df d$ 

COS* i 



+ 



' de d« dt 

COS* t 



+ 



. .didUdu.^^ 

sin t-r: -TT -r: sin* 6 - ^* <^* <** . 

dd dO dO cos t 

di d^i du , 2 A ^^sin*^ 



de de 
cos i 



+ 



3 cos* i sin* t J v} sin* 0 di^ . , ^ 

p — u sm* 0 

cos* to d0 

j d?i (in du a „• a 

cos* t 






cos* 1 



+ 

+ 

+ 

+ 

+ 



'diydu sa ■ A 
— — cos* 0 sin 0 + 

,d0 d0 



cos* t 

4 sin*t^^u cos^0 sin 0 



cos* t 



2 sin i|0§cos ^ sin g 
cos* t 

cos* t 

2sit^*^f cos ^ sing 
cos* I 

d^i dQ du A • A 

2 sin^z^^cos 6 sin 6 sin x-^^cos 0 sin 0 



(i^ 

COS i 



cos X 



+ 



2s'n’‘(f)’^cos«sm« ( di\‘ du . 



cos* t 



di dQ ^2 sin t’4u cos 0 sin 0 

+ 4 sin i— —u cos g sin g : 

d0 d0 cos 1 
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+ 



+ 



+ 



3sinS(§)'ucos^<? 



cos^ i 



cos^ i 



de de ^ _i_ 6 sin cos 0 



cos^ i 



+ 



cos^ t 



, .dt dCl du 2 . 3 sin i^^cos^O 



COS t 



— 3 sin 



in^-^ 



<^Y «d . 2 sin’ igu cos’ « 3 (g) ^ 

— u cos P H ^2—: ^ ^ 

do ) cos t 



cos* i 



+ 



qin£H^ <i^n (fu 
^ d» d» (<»» d» de‘> de 



cos t 

2 _;_ • di dn dti 
Sin t-j7-rTr:i 



cos* i 



cos t 
di dCt du 



cos* t 



cin 

dg dg de _ • ;211 mi — _L ^ d^ de 

olli I , * , - , - l” , 

d^ dO dO cos t 



+ 



cos * iJ u* 
cos* t'o 



+ 



cos* i 
cos* t'o " 



(6.24) 



Equations (5.23) and (5.24) can be greatly simplified if only a first order 
approximation to the equations of motion is desired. The approximate equations 
to (5.23) and (5.24) are respectively: 



-2 sin t— - sin 6 
dO 



dH . .d*n di 

-TTr sin 6 + sin t— — cos 0 — 2— cos 6 
dO^ dO^ do 

2 cos^ t sin i J u sin 0 

cos* t'o 



(5.25) 



^ sin t^g^sin(2g) 

dO^ cos i 



O o I -n 2 • dfl d^u J d^Q du _ • 2 /j _ • * di du _ • 2 >3 



COS t 



cos* t 



cos t 



3 cos* t sin^ t Ju* sin^ 0 



2 sin^t^^cos 6 sin 0 



cos* t 



cos to 

sin t^% cos ^ sin 6 2 sin t^u cos 0 sin 0 



cos t 



cos i 
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(5.26) 



3 sin cos* ^ ^ ^ 



cos t 



cos I 



+ 



+ 



ndn^ <f»n du 



, g.... I COS*tJu* 

cos t cos t 



cos t 



+ 



cos* to 



cos* i 

cos* t'o 



Equations (5.25) and (5.26) will be used in the initial analysis; the general 
equations will be used for second order calculations. 



F. THE INITIAL CONDITIONS 

It is desired that the solution to the equations derived in the last section 
be an accurate, long term predictor of the satellite’s motion. In fact, as long as 
the oblateness perturbation remains the dominant disturbing force, the solution 
should be valid for close to 1000 revolutions. However, before the solution can be 
used for prediction, a set of initial conditions must be determined. 

The subject of orbit determination represents a separate discipline within 
celestial mechanics, and a discussion of the various techniques used to determine 
an orbit is outside the scope of this analysis. It is sufficient therefore, to state that 
the purpose of orbit determination is to find the orbital elements of a satellite from 
reduced observational data. A set of observations will determine the osculating 
elements at time to. As was noted in Chapter 4, these elements will change, and 
at ti a new set of osculating elements may be calculated. For the purpose of this 
analysis any observed set may be designated as the osculating elements at to S’lid 
thus prescribe the initial conditions. Stated mathematically, the initial conditions 
are at t = to: 



a(l — e*) dr _ o(l — e*)c sin(^o “ ^^) 

1 + e cos(^o “ (1 + e cos(0o — <^)* 
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where a = oq, e = eo, io = wo, 



and i — to fl = flo y = ~ ^ 



( 5 . 27 ) 
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VI. THE PERTURBATION PROCEDURE 



A. PRELIMINARIES 

The perturbation method used in this analysis is a variation of the technique 
known zis the Method of Strained Coordinates. The motivation for the use of this 
technique is the subject of this section. 

In perturbation theory, the quantities to be expanded can be functions of one 
or more variables besides the perturbation coordinate. An asymptotic expansion 
of f{0‘,J) in terms of the asymptotic sequence em{J) is, 

f{6;J) ~ f; as J 0 (6.1) 

m=0 

where ^ is a scalar (or vector) variable independent of J and the coefficients 
are functions of 6 only. This expansion is said to be uniformly valid if, 

f{d-,J) = J2<^m{e)€m{J) + Rn{e-,J) 

m=0 

R„{6-J) = 0(£„(J)). 

For these uniformity conditions to hold, a^(^)€TO(J) must be small compared to the 
preceding term a,n-i(^)€m-i(‘I^) for each m. Each term must be a small correction 
to the preceding term regardless of the value of 0. 

Unfortunately, it is the rule rather than the exception that expansions like 
(6.1) are non-\miformly valid and break down in certain cases. A case of particular 
interest is the presence of secular terms such as cos 6 and sin 0 which make 
/m(^)//m-i(^) unbounded as 0 approaches infinity. 

To illustrate how secular terms arise in the solutions to differential equations, 
reference is made to one of the equations of motion derived in the preceding 
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chapter. Equation (5.26) has the form, 



u" + u = 1 + J /(u, u', u", i', i", n', n") . 

When expanded, the right side will be linear combinations of trigonometric terms 
and constants. The presence of terms of the form cos(^), sin((9), cos(3y — 20), 
sin(2y — 6), etc., on the right side will produce secular terms. For example, the 
second order differential equation, 

u" + u = cos 6 



has the solution. 



u = 



0 sin 0 + cos 0 



+ Ki sin 0 + Ki cos 0. 



Note that the presence of the term 0 sin 0 will result in unbounded solutions for 
u as 0 approaches infinity. Since in this representation u is the reciprocal of 
the distance from the center of the planet to the satellite, the physical effect of 
the secular terms would be to produce in r periodic terms with large amplitude 
variations, a situation that certainly does not correspond to physical reality. 

A technique for dealing with terms that produce secular terms is to eliminate 
them from the right side of the differential equation. The Method of Strained 
Coordinates is a perturbation technique designed for removing secular terms. To 
illustrate, it is recalled that the solution to the two body problem is 



l/r = u= l + e cos(0 — w), where p = 1. 

In the Method of Strained Coordinates, the 0 coordinate is strained by introducing 
a new variable y = f [0) = 0 — oj + Jyi0 + .. . As will be shown in the next section, 
choosing the correct value for yj will insure secular terms do not arise in the first 
order solution for u. 
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While the above method removes secular terms to first order, it is not ade- 
quate for dealing with secular terms that arise to second order in the equations of 
motion. That it is important to remove second order secular terms is shown by 
the following equation: 

/(r) = 1 -b e cos(^ — u>) -f J(cos ^ -f . . .) 

-b cos{6 — u;) -b . . .) -b J^[6^ cos(^ - u») -b . 

Note that, although the terms through order J are boimded, zs 6 —* oo, the 
and terms grow without boimd and dominate the right side. In this present 
analysis, 6 has an upper bound of (27t)10^; however, the effect of secular terms 
remains since etc., are all of order J zs 6 (27t)10®. An infinite series 

would have to be retained. 

In this present work, three additional techniques had to be devised to deal 
with secular terms to second order. These techniques will be discussed in the next 
section. The perturbation method therefore is not strictly the Method of Strained 
Coordinates, but a variation of that method. 

The basic steps in the perturbation procedure are as follows: 

1. The dependent variables and independent variables are expanded as func- 
tions of a small parameter (J). 

2. The variables are then substituted into the equations of motion, and the 
equations are solved consecutively. Each solution yields a more exact ex- 
pression for the appropriate variable. 

The process is carried out through second order to demonstrate that all secular 
terms may be eliminated and that the solutions are bounded. The following sec- 
tion highlights the calculations involved in the process and shows the first order 
equations and solutions. The second order expressions are long, and their display 
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in this context would not contribute to the analysis since only a few specific terms 
are relevant. The complete expressions axe contained in Appendix B; they will be 
referred to during the course of the analysis. 

C. SOLVING THE EQUATIONS 

1. First Order Approximation for i and 0 

The equations to be solved are (5.25) and (5.26). Since the right side 
of these equations are analytic functions of J, it is reasonable to expect that the 
solution to u will be arbitrarily close to the two-body solution, 1 + e cos(0 — w), 
when J is sufficiently close to 0. Likewise i, Q, and y will be arbitrarily close to 
to, f^o, and 6 —(x>, respectively, when J is close to 0. This assumption amounts to 
letting 



1 + e cos y + Jui + J^U 2 + . . . 


(5.2) 


0 — u; + Jt/i + J^j/2 + . . . 


(5.3) 


t*0 + + J^l2 + • • • 


(5.4) 


flo + "!"••• 


(5.5) 



The first step in the solution of the equations of motion is to substitute 
(6.2) - (6.5) into equation (5.25) and equate the coefficients of J. The result is 

(— 2^1 sin 6 + flj cos 6) sin t’o ~ 2ij cos 6 — i" sin 6 
= 2 cos t'osin t’osin 0(1 + e cos y) . (5.6) 

Sin t and cos i have been replaced in the above equations by their approximations: 
sin t'o and cos t’o- These are valid approximations since fl". Cl', i" , and i' are all of 
order J. 
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Equation (6.6) is a linear differential equation with two unknown func- 
tions. Terms of the form cos 6, sin 6, etc., on the right side of this equation, and 
the more general equation (5.23), will produce secular terms in t and fl. It should 
be recalled that these same terms ca\ise secular terms in u. Note that (6.6) con- 
tains a sin 6 terms on the right side. There would be a need to eliminate this 
term were it not for the conditions placed on fl and t by the definition of the 
reference plane. Specifically, fl, which governs the rotation of the reference plane, 
must be expressible as an unbounded (secular) term plus bounded periodic terms, 
and t must be bounded. The fact that fl must contain a secular term negates the 
requirement to eliminate the sin 6 term. Therefore, a solution of the form 

fli = ai6 -|- «2 sin y 

I'l = Pi cos y (5.7) 

is assumed. 

By substituting (6.7) into (6.6) and equating coefficients, the following 
solution for i and fl is obtained: 

t’l = — 2/3 e cos t'o sin t’o cos y 

Qi = —6 cos I'o — 4/3 e cos t’osin y (5.8) 

2. Second Order Secular Terms in i and fl 

Equation (6.8) satisfies (6.6) to order J. The next step in the process 
is to substitute (6.2) - (6.5) and (6.8) into (5.26) and solve for ui. However, 
proceeding with (6.8) in its present form will lead to secular terms in second 
order. A brief paragraph will explain why secular terms are anticipated. 

Success in using perturbation methods requires a certain a prtort knowl- 
edge about the nature of the particular problem one is trying to solve. There b a 
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certain aimount of trial and error involved in the process. For instance, the success 
of the Method of Strained Coordinates depends on the knowledge that a secular 
term will arise in the first order solution to (5.26) and that a mechanism (the 
strained coordinate) can be put in place one step ahead in order to eliminate the 
seculair term. The second order secular terms were found by this same trial and 
error process. 

Returning now to equation (6.8), it was discovered in this analysis that 
the problem terms 



(-45s^ + 28s) 
24 



J^e^c sin(2y — 6) 



( 6 . 9 ) 



and 

-J^s^e^c sin(2y — 3^) 

8 

(where s = sin io, c = cos t'o) • 

appear on the right side of (5.23) when it is expressed to order J^. The appearance 
of these terms that give rise to secular terms was unexpected. They were not 
reported by Brenner because the authors assumed a small eccentricity and dropped 
and higher terms. 

There were several failed attempts in dealing with these new terms be- 
fore an effective measure was discovered. First, an investigation was made into 
the effect of retaining the secular terms. 

The secular terms produced by (6.9) are: 



*2 

CI2 



(15s^ - 14s) 
24 

(155^-7) 2 
12 ^ 



■e^c 0 sin(2y — 26) 
c 0 cos(2y — 20) . 



As was demonstrated earlier, when 6 (27r)l0* these terms are of order J and 

must be retained in the first order solution. As the solution progresses, these terms 
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will continue to produce secular terms with coefficients J*6^, . . ,, *. A 

convergent series representation for these terms was not found. 

The remaining alternative was to alter the form of (6.7) to eliminate 
secular terms. If terms with arbitrary coefficients can be found which when added 
to (6.7) produce terms of identical harmonics to the same order as (6.9), the co- 
efficients may be chosen so that all terms with these harmonics are eliminated. 
Essentially, the only terms that may be added to (6.7) are terms which satisfy the 
original differential equation, equation (5.23). Therefore, one may add homoge- 
neous solutions of the differential equation or arbitrary constants. 

Adding an arbitrary constant to (6.7) has the following effect on higher 
order terms. From equations (5.23) and (5.24), it would appear that only the 
derivatives of i and fl enter into the higher order calculations, and therefore arbi- 
trary constants would be eliminated. This is true for Cl, but not for t. To explain, 
it is recalled that the approximations sin i = sin I'o and cos t = cos I’o were valid 
for the first order approximation of t'l and fli. This approximation was valid since 
all terms were of order J. However, in the calculation for ui there is a term 
cos^ i/ cos^i'o that is of order 1 (Eq. (5.26)), therefore a better approximation for 
i is required. Using (6.7) and the Taylor series expansion for cos t (keeping only 
terms of order J) results in 

cos i = cos(i‘o -h f{J)) = cos t'o ~ sin to f{J) 

= cos t'o — sin i’o(*i) 

Adding an arbitrary constant Ka to (6.8) results in 

cos i = cos t’o + sin io(2/3 J e cos t'o sin I'o cos y -|- Ka) 

(A similar expression is required for sin ») 
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Therefore, a constant added to (6.8) will alter the form of subsequent 
terms, but not in the form required to eliminate secular terms in t '2 and 02- It will 
be shown that a constant like Ka will produce a term of the form KaJ sin^ t’o in ui 
eind terms of the form KaJ^e sin* t'ocos t'ocos y in * 2 , and KaJ^e sin* t’ocos t’osin y 
in fl 2 . Although they cannot be used to eliminate secular terms, the constzmts 
will be needed to satisfy the initial condition, therefore it is essential that they do 
not produce irremovable secular terms to higher order. 

The next alternative for the elimination of secular terms is to add a 
solution to the homogeneous equation of (5.23). While it will be shown that this 
technique is successful in eliminating secular terms in U 2 , it did not succeed in 
removing them in t ‘2 and 02- It failed because the homogeneous terms produced 
new secular terms to higher order. Many various combinations of homogeneous 
solutions were added to (6.7), but all attempts along this line only complicated 
the problem. 

An answer to the question of how to eliminate the secular terms was 
suggested in a report by Weisfield [Ref. 31] on polar orbits. When faced with the 
problem of eliminating secular terms from his equation for A*, Weisfield added a 
term like cos(2y — 20). Now to the particular order to which one is working this 
term acts like a constant in the derivative, e.g., let, 

y = 6 + J 6 



then 

^(cos(2y - 20)) = sin(2y - 20)(2 - 2{l + J)) = 0(J). 
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The attempt was then made to apply Weisfield’s technique to this more 
general problem. (6.7) would then become 

*1 = — 2/3e cos t'osin t'ocos y + K\ cos(2y — 26) + 

fli = ^ cos to - 4/3ccos t’osin y + A’ 3 sin( 2 y — 20) + . (6.10) 

where: Ki = —Kic^cs = Ki cos(2w) + 

Kz = —Kzt^c Ki = — ffjsin(2a;) + Ki 

(It is noted that (6.10) still satisfies (6.6) to order J) 

The added haxmonic terms produce identical harmonics as (6.9) to sec- 
ond order and no other new secular terms. Thus, they allow for the elimination 
of these problem terms. 

3. First Order Approximation for u 

With valid expressions for t'l and fli established, the procedure may now 
be continued with the calculation of ui. Substitution of (6.2) - (6.5) and (6.10) 
into (5.26) and letting sin t’o = 5 and cos t’o = c results in the following equation 

d^U\ 5c^s^ cos(2y + 20) 

+ ui = 2ecosyyi-f 



d02 



8 



"I” 2c^K\S^ cos(2y — 20) 



cos(2y — 20) 



8 



^ 5es^ cos(y + 20) lle^s^ cos(2y) ^ 5e^cos(2y) 



- 3e^5^cos(20) s^cos(20) 

— 5es cos y + 4e cos y H : 1- 



— 2e^KiS^ cos(u;) — 2K:S^ — 



4 

17e^s^ 

12 



5s^ 7e^ 



In the above equation, the cos y terms will produce secular terms in Uj. A choice 
of yi = (5s^ — 4)/2 will eliminate that possibility, y becomes 



y — 6 — ui J B ^ 2 ^^ — d^y^B . 



63 



Equation (6.11) becomes 



^ le^£[2y + 26) ^ 2e^i^^5^os(2y - 2d) 



de^ 



8 

cos(2y — 26) 5cs* cos(y + 26) lle^s^ cos(2y) 
''' 8 3 4 

^ 5c^cos(2y) ^ 3 c*5*cos( 20) ^ «^cos(2^) 



I7^2e2 tie* 7#.2 

— 2e*iiri5* cos(2u>) — 2KiS^ ; h — + 1. (6.12) 

12 2 6 

To solve (6.12), a solution of the form 

Ui = ai cos(2y + 26) + Oj cos(2y — 26) + «3 cos(y + 26) + cos(2y) 

ascos(2^) + ao (6.13) 

is assumed. 

Equation (6.13) is substituted into (6.12) and the coefficients of like harmonics are 
equated. The coefficients are: 

—e^s^ 



ai = 


24 


aj = 


+ D 


03 = 


— 5cs^ 


24 




flls^ 5 


04 = 


~r:r ~ 7 


V 12 6 




/e 1\ 


05 = 


“ 7 + ^ 


\4 6/ 


Oo = 


—2Kie^s^ c 




7 2 

+ + 1 




6 



2 2 
€ S 



17e> 5\ , 



64 



Equating (6.13) for Ui results in: 



Ui = 



_ cos(2y + 2^) 

2^2 



24 



+ 2e^K\S^ cos(2y — 20) 



cos(2y — 20) 5es^ cos(y + 20) ^ lle^s^ cos(2y) 



8 



24 



12 



5e2cos(2y) cVcos(21?) s2cos(2<?) 2 , 

r ze KiS cos(zu;j 



j 17eh^ 5s* 7e* 

- 2ifjs* + hi. 

12 2 6 



(6.14) 



4. Second Order Secular Terms in u 

The above expression for Ui is not complete. Looking one step ahead 
it is known (by trial and error) that the following problem terms arise in the 
differential equation for Uj. 



d*U2 



-h U 2 — 



+ 

+ 



f 75e^s* — 200e^s* + 136c^ 

V 240 15(5s* - 4 )^ 

^ 375eV + (- 480e* - 40e)s* + 136e^ 



240 



cos(3y — 20) 

cos(y — 20) 



15(5s* -4), 



f (225c^s^ — 225c*s* + 2e^) cos(2w) 2e*cos(2o;) 

\ 60 75s* -60 

(45e3 - 550e)s'‘ + (36e* + 488e)s* - 56e^^ 



48 



cos(y) 



(6.15) 



(Only the problem terms are displayed. The complete expression is equation (B.5), 
App. B.) 

The harmonics cos(3y — 20), cos(y — 20), and cos y all cause secular 
terms in U 2 , and therefore must be eliminated. In addition there is another more 
troubling problem with (6.15). That is, certain terms have in their denominator 
(5sin*i'o — 4 ), and if t’o = ±sin"^ ^/4/5, then the denominators are zero. This in- 
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clination is well known in Astrodynamics, and has been named the critical inclina- 
tion, or more appropriately the critical inclinations since there are two: t’o = 63°26' 
and 117°34'. The problem of the critical inclination will be dealt with later in the 
analysis. In addition, there is qualitative discussion of the critical inclination in 
Appendix C. 

The task is to eliminate the three terms which give rise to secular terms 
in uj from the right side of (6.15). By inspection it is seen that a proper choice of 
t /2 will eliminate secular terms produced by cos y. It was then discovered in this 
analysis that the addition of a term of the form cos(y — 2^) — cos(y + 2a») eliminated 
the cos(y — 2B) problem term. This term may be added since it is a solution to 
the homogeneous equation. There remained one term to eliminate; cos(3y — 26). 
This term can be eliminated by adding a term with the harmonic sin(2y — 26) to 
yi. In addition, a constant weis added to yi to satisfy the initial conditions. The 
complete first order expression for y is now 

y = Q - (jj j ^ ^£^2 -2^ + sin(2y - 26) + Ke'j 

. (6.16) 



where Ki = — J/fsandXe = — J/i’s sin(2w) -f K&. 

By adding —K 7 (cos(t/ — 26) -fcos(y-f2w)) to eliminate a secular term and K$ cos y 
-fXgsin y to satisfy initial conditions, ui becomes 

cos(2y -f 26) 



Ui = 



24 



+ 2e^Kis* cos(2y — 26) 



cos(2y — 26) 5es^ cos(y + 26) 

-t- ^ ^ - eKr cos(y - 26) 



8 



24 



T, , ^ . lle^s^ cos(2y) 5e^cos(2y) , 

-f eKj cos(y -|- 2u) + — + Kg sin(y) 



12 



T- , , e‘^s^cos(20) cos (2 6) , 

-I- Ks cos(y) 2e^KiS^ cos(2w) 
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(6.17) 



- 2X28"^- 



12 



55^ 





1 . 



5. Second Order Solution for t and fl 

With all terms in place to deal with secular terms the calculations are 
continued by substituting (6.2) - (6.5), (6.10), (6.16), and (6.17) into (5.23) to solve 
for Ct 2 and t’ 2 , and to evaluate the constants Ki and K^. The result is contained 
in Appendix B, i.e.. 



(— 2 fl 2 sin 6 + flj cos 6) sin »o — 2»2 cos 0 — ij sin 6 
= Right side (B.l), (Appendix B). 



An inspection of (B.l) reveals that the coefficients of sin(2y — 6) and sin(2y — 30) 
form respectively the following simultaneous equations; 

{SKs + lOKi - 45)e^s3^ {4K3 + 4Ki -7)ehc 

24 6 “ 

+ (4R'i - 4if3)e^s c = 0. 



Solving these equations results in 



Ki 



Ks 



15s^ - 14 

24 ( 5 s 2 - 4) 

75s* - 120s‘ + 56 

24 ( 55 ^ _ 4)2 



(6.18) 

(6.19) 



Substitution of these values for K\ and Ks into (B.l) gets rid of all sin(y — 20) 
and sin(2y — 3^) terms. With an eissurance that no secular terms will arise, the 
equation for U 2 and (2 can be solved. 

Before progressing, it is noted that (6.18) and (6.19) contain the same 
problem denominator that was observed in (6.15). In fact, (6.18) and (6.19) are 
the first occurrence of the critical inclination term in this analysis. The term then 
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continues to manifest itself in all higher order analyses. The reason for the presence 
of the term is simple to explain. The constants Ki and were multiplied by the 
derivative of j/i “ ^)) <^^ring the analysis. The derivative then shows 

up in the denominator of these constants when they are evaluated. It will be 
important to show later in the results that despite the occurrence of an apparent 
singularity, the final solution is uniformly valid for all inclinations, including the 
critical. 

With secular terms removed the following equation can be solved for flj 

and tj. 

(—202 sin 6 + O2 cos $) sin t'o — 2t2 cos 6 — t" sin 6 (6.20) 

= Right side (B.2), (Appendix B.). 

As was done in (6.6), a solution which includes the harmonics of the right side of 
(B.2) is assumed for fl: and t'2. This solution is substituted into (B.2), and the 
coefficients of the various harmonics are equated. Once again the conditions are 
that fl2 must be expressed as a secular term and bounded periodic terms, while 
t2 must be bounded. The solution is 

fl2 = Right side (B.3) 

t2 = Right side (B.4). (6-21) 

6. Eliminating the Final Secular Terms 

To complete the solution the remaining constants must be found. Three 
of the remaining six constants, 1/2, K5 and Kjy are obtained from the differential 
equation for U2. It is recalled that these constants are used to eliminate secular 
terms in U2. There is no need to solve the equation for U2 since the constants may 
be evaluated from the right side of the differential equation. Therefore, (6.2) - 
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(6.5), (6.10), (6.16), (6.17), and (6.21) are substituted into (5.24). The resulting 
equation is (B.5). 

As in the previous solution for ri2 and t’2, the secular terms in U2 are 
eliminated by setting the coefficients of the problem terms, i.e., (6.17), equal to 
zero and solving for the constants. The results are: 



where 



K, 

K7 

V2 



75e^s^ - 260e^s^ + 296e^s^ - 112e^ 
192 (|s2 - 2)^ 

15c^s^ — (14e^ — 2)s^ 

48 (|s2 - 2) 

e^cos(2w) 



( 6 . 22 ) 



y2 = — lOes^cos(u) — 0o) + 



26cs^ cos(u; — 60 ) 15e^s^ cos(2u>) 



8 



15c^s^ cos(2w) e^cos(2u;) 15e^s^ 275s^ 3c^s^ 61s^ 

+ — — - + ^:z ^ + ^r- + 



8 



60 



32 



48 



8 



12 



If! 

12 



With the above constants evaluated, it is assured that the second order expression 
is free of secular terms. 

7. The Initial Conditions for y, j, d, Cl 

The task now is to evaluate the remaining constants by establishing the 
initial conditions. At t = <0, it is required that the velocity vector of the satellite in 
the reference plane be tangent to the corresponding two body ellipse determined 
by the satellite. In addition to the initial conditions established for r and ^ in 
(5.27), it is recalled that at < = to 



y = 60 — u, i = i‘o, and Cl = CIq 
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From (6.16) and (6.22): 

y = 0 -u + J (e - 2^ 

(75c^ 5® — 260e^5'* + 296e*5^ — 112e^) . \ 

- ^ ; Tj ^(sin(2y - 26) - sin(2w)) 

1Q9(5.c2-9» ^ 



+ Ke- 



192 (§5* - 2 )^ 



(6.23) 



30 (§52 - 2 ) 

Choose Ks such that ai t = to, y = Oq — u; therefore, 

Ke = -J60 (^ 5 * - 2 ) . (6.24) 

To obtain the initial condition for i, use (6.10) and (6.18) so that 



,• ,• -14)(cos(2y-2^) -cos(2w)) 

t — to 2JoJccs cos y Jc cs 24(5 ^ 4) "^2' 



(6.25) 



It is evident that the Je^ terms are eliminated when y = 6 — u> and 6 = 6 q. The 
result 



Ki = 2fZJecs cos (00 — w) 



(6.26) 



gives the desired i = t’o at t = to> 

For fl, (6.10) and (6.19) are required. In addition all secular terms from 
fljj equation (B.3), App. B, are needed since all these terms are of order J. 
Therefore, 



fl — flo — J 6c -(■ 



ce^J{75s* — 102s^ + 56)(sin(2y — 20) + sin(2u>)) 
24 ( 55 ^ - 4)2 



AceJ sin y ce*J* cos(2w)0 ^ 5cc^J*s^ cos(2w)0 



ce^J^ cos(2u>)0 
12 

ce^jH cJ^e 



+ 



15s2 - 12 
+ ZcJ'^KisH - 

+ /^4 



8 

5ce^Jh^6 5cJ^s^6 
24 ^ 3 



(6.27) 
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and 



K4 = JcOq + -Jce sin(0o ~ • 

o 



(6.28) 



8. 6 As a Function of Time. (Initial Conditions for r, 

To complete the evaluation of the remaining constants in u, it is nec- 
essary to give ^ as a function of time. The expression for the perturbed ^ then 
can be related to the well-known two-body formula, and from this relationship, an 
expression for r and ^ may be derived. The task will be to choose Kg and Kg so 
the conditions as given by (5.27) are satisfied, i.e.. 



/ X _ fl(l - e^) ^ _ g(l - e^)e sin(^o - a;) 

° 1 -f e cos(^o — (^) do ^ (1 4- e cos(0q ~ w))* 

Proceeding in this manner, the formula for ^ can be obtained from the 

relation 

dt _ d<l> dt d<f> cos^ 6 

do do d(f> dO h cos to 

Direct substitution of the derived expression for ^ equation (5.21), results in 



dt 

To 



r 

h 



^ , j( c* cos(2y — 20) — cos(2w) ^ cos(2y — 20) 

1 + I . _ I 



15(5s2 - 4) 



8 



+ 



cos(2y — 20) es^ cos(y + 20) cs^ cos(y — 20) 

60 6 2 
4es^cosy 4e cos y s^cos(20) 2es^ cos(w — 0o) 



e* 5 ^cos( 2 o;) c*cos(2u>) 

8 ^ 60 ^ T ” 



(6.29) 



At t = to 



dt 

dO 



(M = 
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The above condition is an integral from the two-body problem. If (6.29) is evalu- 
ated at ^ = 6q, then may be expressed as 

where 

es^ cos (a; + 6q) 2es^ cos(w — 0q) 4c cos(w — 9o) 

= 2 3 3 

cs^ cos(u> — 3^o) s^cos(2^o) 

6 ~2 '^~ 2 ~ 



From this expression, a formula for h results: 



h = /i(l -|- J Kiq) 



so that from (6.29) 



Using 



P = 



dt , , 



'h? h? 



GM GM 



(l -|- 2J i^io) — o(l ~ ^*)(l + 2^ Kio) 



a new formula for r, which includes the first order perturbation effects, may be 
written as 

« a(l - c*) . . 

' ^ (6.30) 



r = - = 



u 1 + e cos y -|- J{—2Kio + Uj) 
where Uj is (6.17). 



The formula for ^ is 



dr _ a(l — c^)c sin y 
d6 (1 e cos y)^ 



+ 



Jg(e^ - 1) 
(1 4- e cos y)2 



sin(2y -h 20) 5es^ sin(y -|- 20) 
6 8 



— eKj sin(y — 20) 
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' ^ 6 3 2 

T^r . , o • V e*s’sin(20) 5’sin(20) 

— K^ sin y + 2e sin y + K<i cos y H — - H ^ — - 

2 3 



(6.31) 



Next, evaluate (6.30) and (6.31) at 0 = ^o- Keep only the terms to order J. The 
result is two simultaneous equations that can be solved for K% and K^: 



K. = 



eV^os(3w - ^o) lle*5* cos(3w — ZBq) 

16 24 

5e* cos(3w — 30o) cos(3w — 5^o) , 31es^cos(2w — 26q) 

12 16 12 

7e cos(2w — 26q) 3es® cos(2w — A6q) lle^s^ cos(w + Oq) 

3 8 16 

cos(w + 0o) lle^s^ cos(u; — 6 q) 7 s^ cos(o; — 6 q) 

H 1 1 r 



8 



31e^cos(o> — 0 q) „ / os 17e^s^ cos(w — 30n) 

— — 3cos(w — 0o) 



12 



48 



7s^ cos(w — 30o) es^cos(2w) 5e5^cos(20o) 13es^ 7e . 

12 2 i + ^ 



sin(3c<; — 5o) lle^s^ sin(3w — 30o) 5e^(3w — 30 q) 

Ks = 7-:;: :r: 1 



+ 



+ 



16 24 12 

sin(3w — SBq) 31es^ sin(2u; — 20 q) 7e sin(2w — 20o) 
16 12 3 

3es^ sin(2ui — 40o) 7e^s^ sin(u; + Bq) sin(w + 0o) 



8 



16 



+ 



67e^s^ sin(w — ^o) 7s^ sin(w — 0 q) , 29e^ sin(ai — 0 q) 

24 2 ^ 12 

, lle^s^sin(w — 3^o) 7s^ sin(w — 30 q) 

+ 3 sin(w - Bo) + + ~ 



48 



12 



+ 



es^sin(2w) 3cs^sin(20o) 



(6.33) 
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Substitution of the values of ifg and Kq in the equation for u insures 



that the remaining initial conditions for r and ^ are satisfied. 

C. THE RESULTS 

1. Preliminaries 

It was noted that throughout the calculations in this chapter the term 
appeared in the denominator of terms in ui, t' 2 , flj, and uj. It would 
appear that the solution is not valid when »o = ±sin“* ^^4/5. And if the solution 
is not valid near these “critical inclinations” , then one should question the validity 
of the underlying perturbation process. 

The purpose of this last section is to display the final solutions for y, *, 
n, and u, and to show that each of the apparent singularities can be eliminated 
in the limit as as fo approaches ±sin~^y^4/5. It is a remarkable eispect of this 
analysis that every term that appears to cause a solution to blow up is exactly 
canceled by a corresponding term that is “hidden” in the solution. 

2. First Order Solution for y 

By use of trigonometric identities, (6.23) may be rewritten as 

y = 6 — u + J6^s^ — 2^ 

-(75s® - 260s“ + 296s^ - 112) / , \\ . / /5 , N\ 

■ o;rL._,V 'OS (2., - je (-. -2))sm(j«(-. -2)) 



96 -2)‘ 

30 (^s2 - 2 )^ ^ - 2 ) + J^0y2 + 0( + . . .) 



(6.34) 



Equation (6.34) is the complete first order solution for y. If the initial 
inclination is such that ||s^ — 2| < 10~®, then (6.34) should be replaced by its 
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limit: 



e*J’sin(2w)0^ SeJ^ cos(u> — 6 q)6 17c^J* cos(2w)0 

y = 0-wH i — 1 i H i — '— 

15 15 60 



+ 



e^J^e 2JH 



60 



+ ^ + o (j’, j‘e, j'e> (!«’ - 2) + . . .) . 



3. First Order Solution for t 

As was done with the expression for y, i is rewritten as: 

i = to + J cs ^^c(cos(0o — w) — cos y] 

c^(15s^ — 14) sin ^2w — J6 — 2^ j sin ^J6 ^|s^ — 2^j 

24 (|s2 - 2) 

+ o{j^ + j^e + j*e^ + ...). 



(6.35) 



Again, as in y, should »o he near the critical inclinations, ||s^— 2| < 10 
then (6.35) should be replaced by: 

4Je, J*0e^sin2w 

I = »o + — [cos(0o - w) - cos(0 - w)J H 



+ O (^sin"to-2) +...) 



4. First Order Solution for Cl 
Rewriting (6.27) for Cl results in 



Cl = Clo — J6c 



+ ce^J 



cos 



(75s^ - 1205^ + 56) 

48(|s2-2)^ 

2w - Je - 2) j sin (^je - 2)) 



4ceJ sin y ce^ cos(2w)5 ^ 5cc^ cos(2w)^ 
3 15s2 - 12 8 
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ce^J^ cos(2u;)0 



12 



+ bcJ'^KiS^e - 



Sce^J^s^e bcJ^sH 



24 



+ 



+ ^ + if, + 0( J’. J’«, J‘«’ + . . .) 

o 2 



(6.36) 



where 



K2 — -Jecscos{0o — u) 
o 

4 

K4 = J c6o + -J ce sin{6o — uj) . 

o 

If ||5^ - 2| < 10-^ then 

r» Tfl 4ce J sin(^ — w) ce^J^sin{2u})6^ Sce^J^ cos{2(jj)6 

n = n „- j«c + + 



+ AcJ^KiO- 



ce^JH llcJ^d ,, 

— — + ^4 



+ o(^j^ + J^e + j*e^ (^^ 5 ^ - 2 ) + . . .) 



where 



K2 = Y5*^^ cos(^o “ 



K4 = \/l/5 JOq + Je sin(0o - w) . 

O 



5. First Order Solution for u 



1 + e cos y 



The complete first order solution for u is: 
e^Js^ cos(2y — 26) 



24 



+ 2e^JK\S^ cos(2y — 26) 



Js^ cos{2y — 26) 5eJs^cos{y + 26) lle^Js^ cos(2y) 

+ 1 eJif, cos(y-2«) 



+ eJifr cos(y + 2a;) 



5e^J cos(2y) e^Js^ cos(20) J5^cos(20) 

6 4 6 

Jifg sin(o.’ — ^o) + JKs cos{u — 6q) — 2e^JKiS^cos{2u>) — 2JK2S^ 

ITcVs^ 5Js^ 7cV 



12 



+ 



+ J + 0{J^) 



(6.37) 
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where 



Kr 

K2 

y 



Kr 

K, 

K, 



15s^ - 14 
48 (I s2 _ 2 ) 

2ec$ cos (00 — w) 

e-u + je - 2) 

^ 2 ( 755 ® - 2605^ + 29652-112) /S , \\ . /_/ 

Je ^cos 2o>- J0 -5^ - 2 ) ) sin ( J0 

48 (§ 52 - 2 ) V \2 JJ \ \ 

e^J^O cos(2w) + J0Q ^ 2 ^^ ~ (see equation 6.22 for yj) 

— 15e25^ — (14e2 — 2)52 

48 (§52 - 2 ) 

Equation (6.32) 

Equation (6.33) 

Applying the condition ||52 — 2| < 10“®, the expression for u is 



^ ^ e^J cos(40 — 2u>) e^J cos(2w) eJ cos(30 — u) 



30 

e^J cos(20 — 2u>) 
10 



"^10 6 
+ J Kg sin(0Q — ^) "1“ K^ cos(0Q — ^) 



/e2j®sin(2u;)02 SeJ^cos(u) — 
+ e cos i — h ^ 



+ 



15 15 

17e2j2cos(2w)0 c2j2 ^ 2J20 



+ 



+ 



6 — bJ 



60 60 5 , 

e(8e2 — 8) J20 sin(0 + w) e^J cos(20) 2J cos(20) 



120 



15 



- +^y.-j + 0 (j\J^e, J®02 [£^2 _ 2 ) + . . .) . 

15 5 30 V \2 J J 

Where K$ and Kg are now; 

e2 cos(3u; — 0 q) e2cos(3u> — 30o) e2 cos(3c<; — 50 q) 



K, = 



20 



20 



20 
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to I cn 



4c cos(2w — 2^o) 3c cos(2w — 4^o) llc^cos(u; + 0 q) 
15 10 20 

cos(ct> + ^o) 89c^cos(w — ^o) cos(u> — ^o) 

5 60 5 

17c^cos(u; — 3^o) 7 cos(w — 3tfo) 2c cos(2a>) 

60 15 5 

- f + 



— c cos(2^o) 



c^sin(3w — ^o) c^sin(3w — 3^o) c^sin(3w — 5^o) 

~ r;: 1 r:: 1 



20 



20 



20 



4c sin(2w — 2^o) , 3c sin(2w — 40 q) 7c^sin(w + ^o) 
+ rr 1 



+ 

+ 



15 10 20 

sin(u> + ^o) llc*sin(a; — tfo) sin(u; — tfo) llc^sin(w — 3tfo) 



+ 



+ 



5 60 5 

7 sin(u> — 3^o) 2c sin(2u>) 3c sin(2tfo) 



+ 



60 



15 



+ 



+ 0(J*) 



The results obtained thus far may be used to predict the orbit of a satel- 
lite for up to 1000 revolutions when a valid set of initial conditions are provided. 
The initial conditions will usually be given as the initial displacement vector r 
and the initial velocity v, whereas the solution obtained is in terms of the orbital 
elements. It is important to show that the coordinates and velocity components 
can be easily recovered from the orbital elements. Indeed, Reference [4] contains 
the criticism that Brenner and Latta did not show how r and v would be ob- 
tained, and that “the derived elements are not such that the velocity components 
are readily obtainable.” (Arsenault et al.. Ref. 4: Vol. 2, p. 5). 

Again referring back to Figure 5.4 and equations (5.1) - (5.3), the po- 
sition of the satellite in the coordinate system may be derived from its direction 
cosines: 



r = r cos 8 cos (j) \ + r cos 8 sin 4>] + r sin 8 k 
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where r is equation (6.30) and zire unit vectors. 



Measuring from the line of nodes results in 

r = r cos S cos(<^ — n)i' + r cos 6 sin(<^ — fl)j' + r sin 6 k 
where i' = cos n i + sin H j 

j' = — sin fli + cos nj. 

Using the angle relationships from equation (5.1) - (5.3), the above expression is 
rewritten as 

T = r(cos^ cos n — cos t sin 6 sin n)i 



+ r(cos 6 sin fl + cos i sin 6 cos n)j + r sin i sin 6 k . (6.38) 



The velocity is found by differentiating (6.38) with respect to $ and noting 



dr dr dO 

di d6 dt 



The result is 




where 



dr I ■ . . n ■ • /> 

— = — cos t sin n— sin 6 — cos tcos n-— r sm 0 

de de de 



+ sin t— sin fl r sin 0 — cos fl r sin 6 
du 



+ cos fl— cos 6 — sin fl— r cos 6 

de de 



cos t sin fl r cos 0 
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.di dO, 

sin t— cos n sin ^ + cos fl-— cos 0 

de de 



+ 

+ 



cos tcos 



n cos 6^ 



dv 

—(cos t cos n sin ^ + sin n cos 6) 
dd 



dr 3 . . - ... .dr . . .dt 

—— = r sm t cos ^ + sin t sin 0— +r cos t sin 0—. 
d0 d$ de 

The requirement now is to find expressions for and the expression for 

^ is equation (6.31). 



Differentiating the expression for i, equation (6.35), and the equation 



for n, equation (6.35), the following equations for ^ and ^ are obtained: 



dt 

d? 

de 



= 2Jecs sin y 



—Jc Jce sin y 

3 ^ 



(6.39) 



The expression for ^ is equation (6.29). By expanding this expression, 
the following equation for ^ is obtained: 



de _ h T ( cos(2y — 29) — e* cos(2o;) cos(2y — 29) 

~di ~ I V 15(5s2 -4) 8 

e^cos(2y — 2^) es^cos(y + 20) es^ cos[y — 29) 

60 6 2 
^ 4es^cosy 4e cos y s^cos(20) 2es^ cos(u; — 0 q) 



c*s^cos(2w) e^cos(2w) s* 

8 ^ 60 "2 ~ V 



For ||5* -2| < 10-3 
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dt 



1 + J 



2e cos(30 — u;) 2e cos{6 + u) 4e cos{6 — u) 2 cos(2^) 



15 



+ 



+ 



15 



+ 



eJ sin (2a;) g 8 e cos(w - ^o) 



15 



+ 



15 



+ 



II 



D. VERIFICATION OF THE RESULTS 
1. Comparison with Brenner and Latta 

The comparison of solutions obtained by Brenner and Latta must be 
restricted to the J, Je, Je^, and terms. Brenner and Latta limited their 
analysis in J to these terms and thus avoided the secular terms that arise in J^e^. 
No mention is made of the critical inclination in their work because ||s^ — 2 | does 
not appear as a divisor in the terms they included. Similarly, this present analysis 
neglected all harmonics in the potential except J 2 while Brenner included some of 
the terms associated with the J 4 (or Di by their notation) harmonic. 

Before comparing solutions, it should be noted that there is some dif- 
ference between the two works in the notation used, this analysis preferring the 
more standard Astrodynamic notation. Brenner used the co-latitude angle (0) as 
a spherical coordinate while the latitude (5) was used here. In addition, Bren- 
ner used M as the independent variable where M + tt/2 defines the angle from 
the ascending node to the satellite. The independent variable 0, measured from 
the ciscending node to the satellite, was used here. Finally, Brenner defined the 
rotation of the line of nodes in an opposite sense to that done here. In summary, 

6 ;r /2 — 0 (where 0 is co — latitude) 

0 M + 7 t /2 (where 0 is the polar angle) 

n -n 
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There will also be a difference in the constant terms. Brenner chose 
initial conditions that would make the analytical solutions as simple as possible. 
This analysis adopted the more general initial conditions that were outlined in 
Chapter 5. 

Now restricting the comparison to the terms that are left, Brenner’s 
solution for i and fl are respectively: 



i = — Jecs cos y+... 

3 

„ 5 , 2 .^ , 4 , . . . . 

Q = JMc -J^Mcs* + -Jec sm y H — — sin(2m) + 

O o 



A check of equations (6.35) and (6.36) shows there is an exact match of terms of 
order J. The final term in Brenner’s solution for fl is a term and can be found 
in the solution for fl 2 , equation (B.4). Next, Brenner’s abbreviated solution for 
Ui is 

5es^ 5^ 5 

cos(y + 2M) H cos(2M) s* + 1 + . . . 

24 ^ ’ Q ^ ’ 2 

These same terms are duplicated in equation (6.17). Finally, Brenner’s solution 

for yi is |s^ — 2 which is exactly the expression obtained in this analysis. 

2. Comparison with an Independent Analysis of Polar and 
Equatorial Orbits 

Appendix A contains an independent analysis of the polar and equatorial 
orbits. Since for the polar case there is no variation of inclination, no rotation of 
the line of nodes, nor dependence of the motion on the longitude 4>, new equations 
must be derived. The equations cannot be solved in terms of the angles i and fl. 

Instead, the equations are solved in terms of the variable which is related to 

it 
de ■ 

For the equatorial case the inclination remains constant. The line of 
nodes is undefined since the satellite remains in the equatorial plane; however, the 
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angle fl does change with time. Instead of using the polar angle 6 to measure the 
angular displacement of the satellite, the angle <f>, measured from the fixed axis 7, 
is used. A new dependent coordinate is defined as 

Y = <f>-K + JYi<l> + J^Yi4 > . 

As Appendix A demonstrates, there is exact agreement between the 
special polar case and the general caise for to = In addition, when the appro- 
priate change of variables is made in the equatorial czise, it agrees exax:tly with 
the general case at t'o = 0°. 
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VII. CONCLUSIONS 
AND RECOMMENDATIONS 

The results achieved in Chapter VI represent a unique solution to the problem 
of a satellite in orbit about an oblate planet. In fact, prior to this present work 
the problem had not been solved in this representation [Ref. 18, p. 340]. In 
contrast to the universal solution presented here, all other methods require a 
reformulation of the problem in alternate variables in order to achieve solutions 
at various singularities, the most prevalent being the critical inclination. 

The results support the theory that at the critical inclination there are no 
discernable features in the perturbations of the elements that distinguish it from 
any other. This conclusion was reached by Lubowe [Ref. 38], is shared by Taff 
[Ref. 18] and has been corroborated by the physical data. However, all analyti- 
cal solutions arrived at via canonic transformations predict perturbations in the 
vicinity of the critical to be 25 times greater than those away from it. 

The perturbation method used in the analysis embodied the principles out- 
lined in Chapter I. First, the resulting solutions are significantly more accurate 
than the two-body solution. For the appropriate orbits, the relative error of two- 
body solution is 1000 times greater than the solution obtained here. Second the 
solution was obtained in parameters closely related to the classical orbital elements 
and in Cartesian coordinates; no transformation to an alternate non-physical set 
of elements was required. Therefore the physical effects are easily distinguishable 
throughout the analysis. Finally, as has been noted the solution is valid for all 
orbital parameters. 
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While the perturbation method was similar to the Method of Strained Co- 
ordinates several novel techniques for dealing with secular terms and apparent 
singularities were introduced. The verification of these techniques was possible 
due to the extensive use of the symbolic computer program MACSYMA which 
allowed for the investigation of higher order formulas. 

The analysis suggests several au'eas for further research. The areas include; 

• The addition of the higher harmonics of the gravitational potential, e.g., Jj, 
J4, etc., to the problem. 

• The investigation into the feasibility of including non-conservative perturb- 
ing forces such as drag or solar radiation pressure. 

• Numerically integrating the equations of motion subject to the initial condi- 
tions used here in order to check the analytic results. 

• An in-depth analysis of the use of canonic transformations in satellite orbit 
prediction to include a comparison with the method used here, the purpose 
being to determine if the critical inclination singularity is an artifact of the 
transformation process. 
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APPENDIX A 



THE POLAR AND EQUATORIAL ORBITS 

A. THE POLAR ORBIT 

Reference is made to the general equations of motion, equations (5.13) - 
(5.15). Now in the case of a polar orbit, the longitude ^ is a constant. Therefore, 
the equations of motion are reduced to 

dh fdeV -GM 3 . J A ,, ,, 

dt^ \dt) r2 V2 2 / ^ ’ 



£ ( 2^ 

dt V dt 



,GMR} . 

— Ji — — — sm(20) 



Let J = ^ u = pjr 



where p = fGM. 



(A.2) 



Define a function A such that 



so that 



2 ^^ Ka 
r — = /lA 

dt 

d hu^A d 



dt 

Applying (A. 3) to (A.2) the independent variable can be changed from t to 6. 
Equation (A.2) may be rewritten as 

dA^ 



dB 



Similarly equation (A.l) becomes 



= — 2Ju sin(20). 



(A.4) 



_ 1 + Ju [t sin(2<?) + 
~dB^'^~ A2 



(A.5) 
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To solve (A. 4) and (A. 5) a solution of the form 



A' = l + JEi + J^Ei + ... 


(A.6) 


u = 1 + e cos y + Jux + + .. . 


(A.7) 


where y = 6 -uj + JOyi + J^ 0 y 2 + ... 


(A .8) 



is assumed. 

Substitution of equations (A.6) and (A. 7) into (A.4) and equating coefficients of 
order J results in: 



dEi 

de 



—e sin(y + 20) + e sin(y — 26) — 2sin(20). 



A solution of (A.9), ignoring terms of order or higher, is 



(A.9) 



_ e cosfy + 20) , 

El = + e cos(y — 20) + cos(2^) 

O 

+ KiC^ cos{2y — 20). 

Substituting (A. 10), (A.6), and (A. 7) into (A.5) yields to order J 



(A.IO) 



dui 

— + ui = 2ecosyyi+ ^ 



- e^Ki cos(2y - 20) 

O 

e^cos(2y — 20) 5ecos(y + 20) e*cos(2y) . . 

+ 1 + : « cos(y) 



+ 



3e^cos(20) ^ cos(20) 1 



The cos y term would produce secular terms in the solution to ux, therefore yx is 
chosen as ^ to make the coefficient of cos y zero. 
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The resulting equation is 



„ _ 5c2cos(2y + 2fi) , c*cos(2y-2d) 

-^ + Ui — g e ill cos(2y — 20j H 

5€cos(y + 2fl) c*cos(2y) 3€*cos(2fl) _ cos(2d) 

+ _ - + _ + _ 

ej__l 
4 2‘ 



(A.ll) 



It is assumed that the solution to (A.ll) has the following form: 

= a© + fli cos(2y + 20) + uj cos(2y — 20) 



+ a 3 COs(y + 20) + cs cos(25). 



(A.12) 



Equation (A.12) is substituted into (A.ll). The unknown coefficients may 
then be solved for by equating the coefficients which have the same harmonics. 
The results are: 



oo 



Cl 



C2 



03 









1 

4 2 

_c^ 

24 

(8Ki - l)e^ 
8 

~24 

12 

1 

4 6 ■ 



Substitution of these coefficients into (A.12) results in 



Ui = - 



€* cos(2y + 20) 



24 



e^Ki cos(2y — 20) + 



cos(2y — 20) 
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^ 5e co s(y + 20) e^cos(2y) e^cos(2fl) cos(2d) 

24 12 4 6 

£!_i 

4 2 ■ 

The following complementary solutions must be added to the particular so- 
lution for u to ensure that any secular terms in the second order solution can be 
eliminated 



—K$ cos(3y — 20) — Kj cos(y — 20) . 

The first order solution for u now has the following form: 

, T /« ,>/l^ cos(2y + 20) 

u = 1 -t- e cos y — J/T6COs(3y — 20) ^ 



2A 



e^J{8Ki - l)cos(2y -20) 5eJ cos(y + 20) 



8 



24 



— J i^7COs(y — 20) 



e^J cos(2y) , ^ 

+ — + J Kgsm y + J Ks cos y 

J. ^ 

(3e= + 2) J cos(2«) eV J , . 

12 r " 2 ■ 

The expressions for Ei and Uj are substituted into (A.6) and (A. 7) respec- 
tively. (A.6) and (A. 7) are then substituted into (A. 4). The result is 



dE2 

d0 



= -/f5sin(3y - 40) -t- 



sin(2y + 40) sin(2y -|- 20) 



24 



+ e^Ki sin(2y — 20) -|- 



sin(2y — 20) 



12 



12 



— e^Ki sin(2y — 40) 



e2sin(2y-40) 5e sin(y + 40) esin(y + 20) 

+ ^ ^ K8Sin(y + 20) + 

-|- Kg cos(y -|- 20) -|- Ki sin(y — 20) H ^ ^ ^ — Kg cos(y — 20) 



— Kj sin(y — 40) -t- sin(3y) -t- e^Ki sin(2y) — 

„ . 5e sin y e^sin(40G'z) sin(40) 
+ Kj sin y r— ^ + r + 



e^sin(2y) 



24 






' sin(20) 



+ sin(20). 
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Integrating this expression for E 2 would yield a secular term in the harmonic 
cos(2y — 20). Therefore the constant Ki is chosen as — ^ in order to eliminate the 
problem term. The equation is now 



dE-, 



FT -inr?n 1/9^ I + c^sin(2y + 2^) 

- -Kb sin(3y - 40) + — — 



dO ' 24 12 

5e^ sin(^2y - 40) 5e sin(v + 4«) _ ^ 

+ ‘ + ifs cos(y + 2«) + Kt sin(y - 2^) + ‘ ~ 

D 2 

- Kq cos(y — 26) - Kj sin(y - 46) + Kb sin(3y) - - 

4 

. 5c sin y e*sin(40) sin(40) 

+ + + -J— 

c^sin(2^) . 

+ |-sin(20). 

Integrating the above expression yields 

P _ V \(f\ e^cos(2y + 4^) c2cos(2y + 20) 

E2 - -if 5 COs( 3 y- 40 ) — + 



48 



+ 

+ 



5e^ cos(2y — 40) e cos(y + 40) _ /Cg sin(y + 20) 

48 24 3 

KsCos{y + 26) ecos(y + 20) ^ 

^ H A9sm(y - 26} 



, , ccos(y-20) i^7Cos(y-4x) ifsCOs(3y) 

+ Kb cos(y - 20) + 



e^cos(2y) 5c cos y c^cosf40) cos(40) 

8 ^ 24 16 24 

c*cos(20) cos(20) 

4 2 ■ 

Continuing the procedure results in the following second order expression 

d}u 2 cKb cos(4y — 26) e* cos(4y — 26) 

d6^ ^ ^ 2 96 

c/fs cos(4y — 40) 5c^ cos(4y — 40) 3c® cos(3y + 40) 

4 ^ 576 16 



(A.14) 
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+ 



+ 

+ 

+ 

+ 

+ 

+ 

+ 



3c*cos(3y + 20) e® cos(3y — 20) 

- 2 X 5 cos(3y - 20) L 

35c^ cos(2y + 40) 5eKg sin(2y + 20) 5eKe cos(2y + 20) 
32 4 + 4 

5e*cos(2y + 20) eif9sin(2y — 20) 

12 4 

eKs cos(2y — 20) eKj cos(2y — 20) eK^ cos(2y — 20) 

4 + 2 + 2 

e* cos(2y — 20) 7c^ cos(2y — 20) 

48 24 

3c/f7 cos(2y — 40) 3cif6 cos(2y — 40) c^cos(2y — 40) 

4 4 32 

c^cos(2y — 40) 7c®cos(y + 40) 7c cos(y + 40) 

4 



8 



8 



5Kg sin(y + 20) _ 5Ks cos(y + 20) 

3 3 

+ cos(y + 20) _ 

16 36 7 yy j 

c® cos(y — 20) e cos(y — 20) 5Kj cos(y — 40) 

16 6 3 

55e® cos(y — 40) 5e cos(y — 40) 5eKs cos(4y) 

144 12 4 

5c^cos(4y) 5i^5COs(3y) 13c®cos(3y) eX9sin(2y) 

192 3 144 2 

ci^gCos(2y) 3 cX 7 cos(2y) 3eK$cos{2y) 

2 4 4 

c^cos(2y) 23e^cos(2y) 25c® cos y 

48 



32 



32 



17c cos y 5cif7cos(40) ^ 5c'*cos(40) 65c®cos(40) 

24 4 192 32 

5 cos (40) 3ci^8cos(20) c/^7cos(20) 

8 2 2 

c^cos(20) 37c^cos(20) cos(20) eKs 

96 48 3 2 

eK 7 5e* 167c® 1 

~4~ "^ 576 288 6’ 



(A.15) 
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It is noted that three harmonics in (A. 15) would produce secular terms in 
the solution to U 2 . The harmonics are cos(3y — 26), cos(y — 26), and cos y. These 
harmonics may be eliminated by choosing 

-3 



K, = -- 

Kr = 

V2 = 



96 

(3c’ - 8c) 
96 

25c’ + 34 



96 



The first order solution for (A.4) is 



a’ = 1 + J (cos (2^) + c cos(y — 26) 

c cos(y + 26) c’ cos(2y — 26) 



12 



(A.16) 



u may now be rewritten as: 

c’J cos(2y + 26) cos(2y — 26) 5cJ cos(y + 26) 

u = 1 + c cos y H ^ — 

^ 24 24 24 

c’J cos(y — 26) eJ cos(y — 26) ^ c’J cos(2y) c’J cos(20) 

+ 24 12 ^ 12 4 



J cos(20) r • / /> \ T t tx \ ^ 

^ — - — J sm(w — ^o) + cos(a; — 60) — 

6 4 2 



(A.17) 



where 



. J6 J s\n(2y — 26) 

y = 6-u + — + 

^ 2 48 



(25el+34) 

96 



Equations (A.16) and (A.17) satisfy the equations of motion to first order, 
and the equations produce no secular terms to second order. 

A comparison may now be made between the polar solution derived in this 
appendix and the general solutions for y and u, equations (6.22) and (6.37) respec- 
tively. To make the comparison easier, the relatively complex initial conditions 
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(equation 5.27) are discarded. Instead, following Brenner and Latta, the initial 
conditions are selected in order to make the expressions as simple as possible. 
In addition, constant terms like sin(2u;) and cos(2u>) are dropped. It will be re- 
membered that these constants were added to the general solution to prevent 
singularities at the critical inclination. It is a simple exercise to add the cos(2tj) 
term to the polar case, specifically to equation A. 10, and then show that identical 
results are achieved with the general case at ig = 90°. In fact, this has been done. 
However, to display these terms adds nothing to the comparison and only serves 
to make the analysis more tedious. 

Checking equations (6.35) and (6.36), it is seen for the polar case, cos ig = 
c = 0, there results 

i = t’o and fl = flo . 



This agrees with orbital theory that for the polar case there is no variation of the 
inclination nor does the line of nodes rotate. 

Setting sin t „ = 5 = 1 and disregarding the appropriate constants, equations 
(6.37) and (6.22) are 

J cos(2y 20) 5e^J cos(2y — 20) 5eJ cos(y -|- 20) 

1 -I- e cos y H ^ ^ ^ ^ 



24 



24 



24 



^ e^J cos(y — 20) eJ cos(y — 20) ^ e^J cos(2y) e^J cos(2^) 



24 

J cos(20) 



12 



12 



— J Kq sin(w — Oq) -|- J K% cos(w — Oo ) — 



eV J 



4 2 



and 



, JO e^sin(2y-20) »J25e^ + 34) 



These results agree exactly with the (A. 17) and (A. 18). 



93 



B. THE EQUATORIAL ORBIT 



Reference is again made to the general equations of motion, (5.13) - (5.15). 
An orbit that remains precisely in the equatorial plane is possible if only even 
harmonics of the potential are considered. Now 6 = 0®, ^ = 0, = 0, and the 

equations reduce to 




(A.18) 

(A.19) 



As in the general case and the polar case the independent variable is changed 
by use of the integral obtained from (A.19) 




£ 

dt 



constant = h 



r d 



Substitution of (A. 20) into (A.19), results in 

d}u 2 

— +u=l+Ju 



(A.20) 



(A.21) 



where : J 

u 



3 

2 




p/r. 



To solve (A.21), assume a solution of the form 

u = 1 + c cos y + Jui + J^U 2 -f . . . (A.22) 

where Y = 0 - u; + + J Vi^2- 



Substituting (A.22) into (A.21) and equating coefficients of order J results in 
d}u\ 



d6^ 



+ ui = 1 + 2J Yie cos Y -f 2Je cos Y + c cos Y. 



(A.23) 
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The presence of the term cos Y will cause secular terms in the solution for ui, 
therefore choose Vj = -1. Equation (A.23) is now 



cos^ 



d^Ui ^ ^ 

d<f>^ 2 2 

Assume a particular solution to the above equation 



(A.24) 



= cto + tti COS 2F. (A. 25) 

Substitution of (A.25) into (A.24) results in the following solution for ui 

1 

ui = 1 + — — -e*cos(2y) . (A. 26) 

2 6 

The approximation for u is: 

u = 1 + e cos y + J ^1 + ^ — ge^cos (2y)^ . (■'^•27) 

The solution must be checked to ensure it causes no secular terms to second order. 
Substituting (A. 26) and (A. 22) into (A. 21) results in 



d^uz 

~d^ 



+ Uz = 



e»cos(2y) 



+ cos(2y) + 2e cos Y Yz 



Se^COSy „ -r 7 r. 

+ h 3 e COS y + H- 2 . 



Note that the terms with cos Y will cause secular terms in uz- Therefore, Yz is 
chosen as: 

5e^ 3 

~~L2 ~ 2‘ 



The first order solution for u is equation (A. 27) 



where Y = 4> — u> — J<f> — J^<f> 



5e^ 



+ x . 



(A.28) 



12 2 ^ 

Equation (A. 27) is the correct solution in terms of the variable 4>. To compare 
the solution found here to the general solution for the case ^ = 0° will require that 
(A.27) and (A.28) be modified such that they are in terms of 6. 



95 



Referring baw:k to Figure (5.1), it is seen that <f> is measured from the fixed 
axis T such that 

<f> = ir/2-\-0 + n. (A.29) 

Now, from equation (6.27), the first order solution for fl with to = 0® (sin to = 
5 = 0 2 ind cos io = c = l) is 

7Jg2 A J^6 

n = n<, - sin(2y - 20) — -Je sin y 1- — — . (A.30) 

48 3 6 2 

(Note: As in the polar case, the constant terms cos(2w) and sin(2oi;) have been 

dropped in the general case for comparison purposes. As explained previously, 

eliminating these terms does not affect the validity of the comparison.) 

Substitution of (A.30) and (A.29) into (A. 28) results in 

7Jg2 7 A 

Y = 6 — u — 2J0 H sin(2y — 20) J^e^0 Je sin y. (A. 31) 

48 ' 12 3 ' ^ 

Now, from (6.22), the solution for y with t’o = 0° is 

7 7 

y = 0 — u — 2J6 J^e^0 H Je^ sin(2y — 2x). (A. 32) 

12 48 

Using equation (A. 32), equation (A. 31) may be written as 

Y = y - ^Je sin y. 



And therefore by use of the Taylor series expansion 



cos(y) = cos y H — Je sin^y. 
3 



Equation (A. 27) may now be written in terms of y{0)'. 



u = 1 + e cos y + J 






2 ^2 

+ -e^ cos 

6 6 




(A.33) 



Reference now is made to equation (6.37), the general solution for u. By 
letting sin t<, = s = 0 in this equation, it is easily seen that the general expression 
becomes (A.33). 
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APPENDIX B 



2nd ORDER ANALYSIS 



Substitution of (6.2) - (6.5), (6.10), (6.16), and (6.17) into (5.23) results in 



. .dCl . . 






dii 



— 2 sin I— sin 6 — sin 6 + sin i cos 6 — 2 — - cos 6 



de 

= ce’XiS*sin(3y — ^) — cK^s sin(3y — 6) + 



de^ de 

ce^Ki s sin(3y — 6) 



Sce^Kis^ sin(3y — se) 



+ cKss sin(3y — Se) + 



cc*ifi5sin(3y — 3^) 



7c e^s^ sin(2y + 3^) 143ce^5® sin(2y + S) 29ce^s sin(2y + S) 

24 ^ 72 18 



+ Sce^KzS^ sin(2y — ^) + lOce^ii'iS* sin(2y — 6) — 



15ce^5^sin(2y — S) 



— 4ce^K3Ssin(2y — S) — 4ce^KiSsin(2i/ — e) + 



7ce^s sin(2y — e) 



— Sce^K^s^ sin(2y — 30) + 

— 4ce^/Ti5 sin(2y — 30) — 



5c€^s^ sin(2y — 30) 



+ 4ce^i^35 sin(2y — 3^) 



35ces^sin(y + 3^) 



24 



+ ce Kys sin(y + ^ + 2u;) 



— 2ce^K\S^ cos(2w) sin(y + 0) — ce^Kis cos(2u;) sin(y + 0) 

5ces^ sin(y + 0) 



— 2 ceK 2 S^ sin(y 4- 0) 



+ 2cKgs sin(y + 0) 



, ... 4ccssin(y + ^) 

— cc/f25 sin(y + ^) ^ — 2cifg5cos(y + ^) 



— ce/f 7 S sin(y — ^ + 2u;) + 



lOcc^/fiS^ cos(2u;) sin(y — 0) 



ce^KiScos{2u>) sin{y — 0) _ 10cci<’25* sin(y — 5) _ 5ce®/fiS*sin(y — ^) 
3 3 3 



41ces^sin(y — 0) 



12 



— 2cKsss\n{y — 0) — ceKrS sin(y — 0) 



ceK 2 Ss'm{y — 0) ce^KiSsin(y — 0) ^ 8cessin(y—0) 
3 6 ’’’ 3 
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+ 2cKgs cos(y — 9) — ce^KiS^ sin{y — 3^) + 

4 



+ ceK 7 S sin{t/ — 3^) — 



ce^KiSsm{y — Z6) 5cc^s* sin(3^) 5c5*sin(30) 



— lOce^Kis^ cos(2u>) sin 6 — 10 cK 2 S^ sin 9 + 

lOcs^sin^ cc^sin^ . . 

- H r cs sin 9. 



Sce^s ^ sin 9 
12 



3 • 3 <B') 

Choose Ki = 24 ^( 6 ^ri 4 ) and Kz = ^^;4^s^3°*4)t^^ to eliminate the secular terms 
sin(2y — 9) 2 ind sin(2y — 3^) in equation (B.l), so now 

, • yi d %2 • y» • *d O 2 yi ^dtz ft 

— 2 sm I-—- sin 9 — sin 9 + sin x cos 9 — 2—— cos 9 

d9 d9^ d9^ d9 

SQce^s sin(3y — 9) ce^s^ sin(3y — 0) , 

= --1155731^ + ^^-cK,ss,.(3y-S) 

llce^s sin(3y — 9) ^ 35ce^s sin(3y — 3^) 5ce^s^ sin(3y — 39) 

240 ^ 1800s2 - 1440 24 



7ce^s sin(3y — 30) 7ce^s^ sin(2y + 30) 
144 24 

143ce^s® sin(2y + 0) 29ce*s sin(2y + 0) 35ccs® sin{y + 30) 



+ cKzs sin(3y — 30) + 

+ 



72 



+ ceK 7 S sin(y + 0 + 2u>) + 



18 24 

78ce^s cos(2cj) sin(y + 0) 



1800s* - 1440 
ce^s^ cos(2u;) sin(y + 0) llce^s cos(2u;) sin(y + 0) 

4 120 

— 2 ceK 2 S^ sin(y + 0) + ^ces sin(y + 0) ^ 2cKgs sin(y + 0) 

8 

— ceK 2 S sin(y + 0) — 2cKgS cos(y + 0) 



3 

70ce^s cos(2u>) sin(y — 0) 
1800s* - 1440 
Sce^s^ cos(2cj) sin(y — 0) 7ce^s cos(2u>) sin(y — 0) 



— ceKjs sin(y — 0 + 2(jj) — 

+ 



12 



72 



35ce®s sin(y ~ ^) _^ 10ce/T2S® sin(y ~ ^) _^ 5ce*s®sin(y — 0) 



1800s* - 1440 



24 
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41ce5^sin(y — 0) 

2cK%s sin(y — 0) — ceKjs sin(t/ — 6) 



_ ceKjS sin(t/ — 0) 7ce^s sin(y — 0) 8ces sin(p — 0) 



144 



+ 



+ 2cK^s cos(y — 0) 



Z9ce^s sin (y - 3^) ceV sin(y - 3^) 9ces^ sin(y - Z0) 

18005^ - 1440 8 4 

+ ceK,s sin(y - Z0] - ^ ~ 

^ ^ 240 4 3 

240cc* 5 cos(2c«;) sin ^ Sce^s* cos (2a;) sin 0 ce^s cos(2o;) sin 0 

+ 1800.» - 1440 + W 



— lOciCjS* sin 0 + 



Sce^s* sin 0 lOcs^sin 0 ce^s sin 0 

- + 



— cs sin 0 . 



(B.2) 



12 3 ■ 3 

Solve for tj and Clj using the technique described in Chapter 6, Section 5. The 
second order expressions for *2 and are respectively 



*2 = 



+ 



+ 



+ 



llce^s cos(3y — 20) ce^s^ cos(3y — 20) 2Zce^s cos(3y — 20) 
900s2 - 720 6 360 

2cKi$ cos(3y — 20) 25ce^s^ cos(2y + 20) 29ce^s cos(2y + 20) 

3 96 144 

llce^s cos(y — 20) 2cKgs sin y 22ce^s cos(2o;) cos y 
900s2 - 720 3 900s^ - 720 

ce^s^ cos(2o;) cos y 23ce®5 cos(2o;) cos y ^ See cos y 
3 180 3 

ce^s^ cos y 85ces^ cos y 2cK^s cos y 2cKjs cos y 
6 36 3 3 

2ceK2S cos y 23ce®s cos y 20ces cos y 



360 



+ 



(B.3) 



n. = 



+ 



+ 



2ce^sin(3y — 2^) sin(3y — 2tf) ^ ce^ sin(3y — 20) 

75s2 - 60 4 30 

4ci^5 sin(3y — 20) I7ce^s^ sin(2y + 20) 29ce^ sin(2y + 20) 

3 72 144 

7ces^ sin(y + 20) 2ce^sin(y — 20) ce^s^ sin(y — 20) 



36 



75s2 - 60 



+ 



12 



3ces^ sin(y — 20) 2 c/f 7 sin(y — 20) llce^sin(y — 20) 
2 3 360 
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4ce^cos(2w) sin y ^ cc^s^ cos(2o;) sin y ce^cos(2w) sin y 
755* - 60 2 15 

, „ 2 . ce^s^siny lOTces^siny 40/^8 sin y 

+ iccK,s s.n y + ^ ^ 

2cKj sin y 23ce* sin y 28cc sin y 4cKg cos y 



360 



+ 



9 



+ 



5ce^s^ sin(20) 5cs*sin(2^) Sce^ cos{2u))d Scc^s^ cos(2w)^ 

16 12 75s2 - 60 8 



ce^cos(2cj)^ 



12 



K r/- 2/1 5ce^s^6 5cs^6 ce^O c6 

+ 5cK2S^e h 



24 



6 



(B.4) 



Substitution of (6.2) - (6.5), (6.10), (6.16), (6.17), (B.3), and (B.4) into equation 
(5.24) results in 

d^U 2 120e^ cos(4y — 20) 5e*s* cos(4y — 26) 

de^ 36000s* - 28800 

leKis^ cos(4y — 20) 31c^s^ cos(4y — 20) 



16 



+ 



+ 



96 



— 3ei^5 cos(4y — 26) 



e^cos(4y — 2^) ^ e^cos(4y — 40) ^ 9e^ cos(4y — 40) 



+ 



+ 



240 ' 15000s< - 24000s* + 9600 36000s* - 28800 

29e^s^ cos(4y — 40) SeK^s^ cos(4y — 40) 923e^s* cos(4y — 40) 

192 4 5760 

cifs cos(4y — 40) 89e^ cos(4y — 40) 3c*s^ cos(3y + 40) 



7200 



16 



e®s^ cos(3y + 20) 21e^s* cos(3y 4- 20) 29e^ cos(3y + 20) 

48 8 12 

480e® cos(3y — 20) 5e*s^ cos(3y — 20) 



36000s* - 28800 
5e*s* cos(3y — 20) 



+ 



16 



+ SKs cos(3y - 20) + 



- 10ii'5S^cos(3y - 20) 

17c* cos (3y — 20) 



30 



35e*s'* cos(2y + 40) ScKts^ cos(2y + 20 + 2w) 

32 4 

3eK^s^ cos(2y + 20 — 2u>) 360e^ cos(2w) cos(2y + 20) 

4 36000s* - 28800 

5e^s^ cos(2w) cos(2y + 20) ^ 19e^s* cos(2u;) cos(2y + 20) 



16 



96 



100 



+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 



e* cos(2u>) cos(2y + 20} cos(2y + 26) QSe^s"* cos(2y + 26) 

80 2 ^ 12 
bt^K^s^ c os(2y + 26) 6be^s^ cos(2y + 26) SeKjS^ cos(2y — 26 + 2u) 
4 8 ^ 4 

3eK$s^ cos(2y — 26— 2u») 4e* cos(2u;) cos(2y — 26) 

4 150005^ - 240005^ + 9600 

172e^ cos(2w) cos(2y — 26) e*s* cos(2c<;) cos(2y — 26) 

360005* - 28800 8 

7c^5^ cos(2o>) cos(2y — 26) ^ 23c'* cos(2o;) cos(2y — 26) 

80 3600 

blbe^Ki cos(2y — 26) 164c'* cos(2y — 26) 2400c* cos(2y — 26) 

360005* - 28800 360005* - 28800 360005* - 28800 

iTj- A e*5* cos(2y — 2^) 257e*5* cos(2y — 20) 

2 72 

5c/^ 75* cos(2y — 20) IbeK^s^ cos{2y — 26) I7c*/i'2 5* cos(2y — 20) 

2 ^ 6 30 

29c* 5* cos(2y — 20) 245c*5* cos(2y — 20) , .. 

60 72 1 \ y ) 

beKs cos(2y — 20) cos(2y — 20) 67e* cos(2y — 20) 

3 50 ^ 3600 

c^cos(2y-20) Oc* cos(2y - 40) 528c^ cos(2y - 40) 

12 360005* - 28800 360005* - 28800 

c*5* cos(2y — 40) 5c*5* cos(2y — 40) SeKjS^ cos(2y — 40) 

2 ^ 16 4 

9cXs 5^ cos(2y — 40) 121e'*5^ cos(2y — 40) 23c^5^ cos(2y — 40) 

4 240 240 

3eifscos(2y — 40) 29c*cos(2y — 40) 11c* cos(2y — 40) 

2 ^ 800 600 

ZeK-jS^ cos{2y + 2u) _ . 3c/rs5* cos(2y — 2w) 

2 2 

, . 7c* 5* cos(y 4- 40) 7c5* cos(y + 40) 

c/fsCOs(2y - 2u) 

7^75* cos(y + 20 + 2u) 7KsS^ cos(y + 20 - 2uj) 

3 3 



101 



2240e® cos(2a;) cos(y + 2d) 5e^s* cos(2w) cos(y + 2ff) 

3600052 _ 28800 4 

7e®5^ cos(2a;) cos(y + 26) ^ 7c* cos(2w) cos(y + 26) 

12 90 



+ 

+ 



4 ! 661c* 5^ cos(y + 20) 47es^ cosfy + 20) 

106^25“ cos y + 20 + ^ ^ 1 

144 36 

lOeifa-s’ cos(y + 20) 131c*5* cos(y + 20) „ i 

^ ^ ^ - 5c 5* cos(y + 20) 

cos(y + 20) ^ 2Kjs^ cos(y — 20 + 2w) + 2K$s^ cos(y — 20 — 2w) 
36 

480e»cos(i,-2<) 25e».<cos(y - 2») , _ 

3600052-28800 7 vv ; 

2e*5* cos(y — 20) — 



16 

3„2 o/l^ cos(y - 20) ^ 17e*cos(y - 20) 



— 8/f7COs(y — 20) + 



1120c* cos (y — 40) 5c*5^ cos(y — 40) 535^ cos(y — 40) 

3600052 _ 28800 ^ 48 12 



+ 

+ 

+ 

+ 

+ 

+ 



5Kjs^ cos(y — 46) 7e^s^ cos(y — 46) 7c*cos(y — 40) 570c^cos(4y) 

3 24 180 3600052 _ 28800 

5c^5^cos(4y) _ SeKsS^ cos(4y) 19c^5* cos(4y) 5cXscos(4y) 

32 4 ^ 192 2 

13c'*cos(4y) 1120c* cos(3y) 269c*5^ cos(3y) 5/fs5* cos(3y) 

80 3600052 _ 28800 48 3 

679c* 52 cos(3y) 677c*cos(3y) _ 816c^ cos(2w) cos(2y) 

72 lio 3600052 _ 28800 

llc^5'^ cos(2o)) cos(2y) 359 c ^52 cos(2cj) cos(2y) 

8 240 

17c^ cos(2u;) cos(2y) 6c'*cos(2y) 528c2cos(2y) 

600 3600052 _ 28800 ~ 3600052 - 28800 

2r^ A \ e‘‘5^cos(2y) 695cVcos(2y) leKjS^ cos{2y) 

llc^/fj5Vos(2y) + — 

5eK^s^ cos{2y) 21c2/Cj5* cos(2y) 121 c '‘52 cos(2y) 

4 2 240 

3007c*5* cos(2y) , , eKsCos(2y) 29c^cos(2y) 

— + eKjCos{2y) + + — 
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1137e^cos(2y) 960e*cos(2u;) cos y ISe’s'* cos(2w) cos y 

200 36000s* - 28800 4 



15e®s* cos(2w) cos y _ e*cos(2w) cos y 
4 30 

15e^s^ cos y 275es'* cos y 



16 



+ 



24 



— 26 eK 2 S* cos y — 



+ 306X25^ cos y 

3e^s* cos y 



+ 

+ 



61es^cos y 7e^cos y SeKjS^ cos{26 + 2oj) eiTss* cos(20 + 2w) 

6 6 4 4 

eKjs^ cos{26 - 2oj) ^ SeK^s^ cos{26 — 2uj) 18Oc^cos(40) 

4 4 36000s* - 28800 

5e^s^cos(40) 65c*s^ cos(40) 5s^cos(40) Seif 7 S^ cos (40) 

— 



32 



8 



+ 



19e'‘s* cos(40) e* cos(40) 432e^ cos(2w) cos(20) 

192 160 36000s* - 28800 

1056c^ cos(2w) cos(20) 3e^s^ cos(2w) cos(20) 



36000s* - 28800 



8 



+ 



e*s^ cos(2w) cos(20) ^ 19cS^ cos(2w) cos(20) 

2 80 

23e^s* cos(2a;) cos(20) ^ 3e'* cos(2tj) cos(20) ^ lie* cos(2a;) cos(20) 



120 



200 



300 



+ 



cos{26) - 4 K, 3 > cos(2«) - l£!£l£2!M 
36000s* - 28800 ^ ^ ^ ^ 16 

127e*s^ cos(20) 7s'‘cos(20) ^ Oeifjs* cos (20) ^ 3e*if25^ cos(20) 

■ H z I" 



16 



+ K 2 S ^ cos(20) + 



6 ' 2 ■ 2 
193e^s* cos(20) 95e*s* cos(20) 2s* cos(20) 



480 



r. 19e^cos(20) 

— 4eKjCos{26) H :;::r7 + 



12 

2e^cos(2w)* 



50e^ cos(2w)' 



+ 



300 ISOOOs^ - 24000s* + 9600 36000s* - 28800 

e‘‘s'* cos(2u;)* 23e^s* cos(2w)* 7e^cos(2u;)* 576e*if2 cos(2w) 



32 



960 



3600 



36000s* - 28800 



304e^cos(2w) 2400e* cos(2w) ^ K^s* cos{2uj) ^ 17e'‘s‘‘ cos(2w) 

36000s* - 28800 ~ 36000s* - 28800 ^ 24 



+ 



5e*s^ cos(2a;) 3eif7S* cos(2w) 3eif6S* cos(2w) 
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19e*/f25*cos(2cL)) 533c^s^ cos(2u)) 41c*s* cos(2ti») 



60 



720 



24 



+ eKjcos{2u) 



A €*^^2 008(20;) 19c^cos(2w) c*cos(2o;) 

+ </fsC<«(2u.) gj + ^i55 



+ 



17e* 



+ 2Kls* + + 20KiS* 



150005^ - 240005* + 9600 360005* - 28800 ' ■ 3 

25c^5^ 437c*5^ , 975^ 7eKjs^ j ZU^K2S^ 

^2 



+ — m 1 ;; 77; — llii' 25 ^ — 



192 



288 



8 



12 



139c*5* _ 2 eKj 17 e* lie* ^ 



847e^5* 

5760 

(B.5) 
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APPENDIX C 



THE CRITICAL INCLINATION 

This appendix contains a discussion of the critical inclination to include its 
mathematical and physical basis and the methods used by various authors to deal 
with the problem. 

In Chapter 6 it was found that the divisor often appeared in sec- 

ond order terms. A brief explanation was given for the appearance of this divisor, 
and it was shown that the apparent singularity may be dealt with in this analy- 
sis in a straight-forward fashion. However, in other methods the problem divisor 
hcis caused enormous complications. This divisor or “critical inclination” problem 
appears to be universal. Hughes’ investigations [Ref. 32] reveal that every ana- 
lytical theory on orbit perturbation analysis contains the critical inclination as an 
inherent characteristic. As stated in Chapter 2, the divisor causes an irremovable 
singularity in Kozai’s method. Its presence has prompted many authors to devise 
very elaborate techniques to get around the problem with the result that many 
otherwise elegant methods are somewhat diminished. 

The mathematical reason for the appearance of the problem term is as fol- 
lows. It is recalled that one of the secular effects of the Earth’s equatorial bulge is 
to cause the line of apsides to rotate. A first order approximation for this rotation 
is given by the rate of change of the argument of perigee which from Allan [Ref. 
33] is 

d'jj _ ZJ-iU 
dt 4{l — 

It is readily seen that if t = 63°26', the perigee position is stationary and the 
line of apsides does not rotate. When the usual Poisson method of successive 
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approximations in terms of a small parameter is adopted for perturbation analysis, 
the above equation appears in the denominator of higher order terms. 

Coffey [Ref. 34] gives an excellent background discussion on the history of 
the critical inclination problem. He states that A. A. Orlov (1957) was perhaps 
the first to notice the unusual situation at the critical inclination. Orlov found 
that the conventional Lindstedt- Poincare algorithm failed at this critical inclina- 
tion and at zero inclination. Nevertheless, Orlov ignored these exceptional cases 
and developed his theories without addressing the problem of singularities. Krause 
(1952), Roberson (1957), and Herget and Musen (1958) had all overlooked this 
difficulty in their assessment of the long term effects of the term on the Kep- 
lerian ellipse before Brouwer (1958) finally pointed out the problem to the latter 
authors. Brouwer hoped that his use of von Ziepel’s method of eliminating vari- 
ables through canonic transformations would allow him to dispose of the terms 
that lead to singularities. However, as was pointed out in Chapter 2, his solution 
contained singularities at, not only the critical inclination, but also zero inclination 
or eccentricity. 

Finding no method to fit the critical inclination problem into current theories, 
astronomers have devised alternative theories to deal with orbital inclinations near 
63°. 4. The critical inclination, they reasoned, was a small divisor problem and 
therefore could be handled in the standard way introduced by Bohlen [Ref. 35). 
That is, in the vicinity of the critical inclination, the solution may be expressed 
in terms of the square root of a small parameter, i.e., y/j instead of J. [Ref. 36). 
The resulting solutions are totally different in character to the non-singular case 
[Ref. 30]. For example, instead of processing secularly the argument of perigee 
will tend to oscillate about a central position. Taff notes the lack of mathematical 
rigor in this method of handling the critical inclination by stating that expansions 
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in y/j and theories about the libration of perigee are a result of misapplication of 
perturbation theories [Ref. 18, p 340). 

Usually when singularities arise in celestial mechanics there is a physical 
reason for their occurrence. For example, singularities in the solutions for the 
motion of asteroids and natural satellites are caused by the faet that their orbital 
periods are nearly commensurate with that of the disturbing body. Under the same 
circumstances, an artificial satellite may encounter resonant perturbations in its 
orbital elements caused by tesseral harmonics. However, such a physical reason 
does not seem to suggest itself for the critical inclination problem. Message [Ref. 
37] implied that there may be a resonance between the satellite’s mean motion 
relative to perigee and its mean motion relative to the ascending node; however, 
most skeptics have not been convinced [Ref. 32). If there were resonance one would 
expect to be able to experimentally measure deviations caused by resonances, 
but satellite guidance engineers have reported no adverse effects at the critical 
inclination. 

Another check on the physical effects of resonance should be a numerical 
check of the equations of motion at the critical inclination; however, this too 
has been inconclusive. Lubowe [Ref. 38] carried out a study in which he inte- 
grated a satellite’s equations of motion for periods up to 24 hours with various 
initial conditions at inclinations near and away from the critical. He found no 
discernible features in the perturbations which distinguished the critical from any 
other. However, in a separate analysis, Hughes came to the opposite conclusion. 
Hughes developed the Hamiltonian in terms of the Hill variables and showed that 
the critical inclination remained an inherent part of the solution throughout a 
multitude of various canonic transformations. He then integrated the equations 
of motion numerically and was able to show some resonant effects in the perigee 
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height for satellites at the critic zd inclination. He then reasoned from this analysis 
that the critical inclination is unique, that it represents actual physical resonance 
and that it is not merely a by-product of the method of solution or the type of 
variable used in the analysis. 

This disagreement characterizes the lingering debate over the question of 
whether the critical inclination is an artifact of the analysis process or a physical 
reality. Most textbooks choose to mention the problem; however, few voice an 
opinion. Roy [Ref. 20] devotes only two sentences to the subject, merely noting 
that some perturbation analyses break down there. Hagihara [Ref. 39], mentions 
that the critical inclination is an essential feature of the perturbation process 
but notes that heated discussions continue concerning its physical reality. Other 
authors, namely Herrick [Ref. 40], Kovalevsky [Ref. 41], and Geyling [Ref. l] 
mention the problem and the standard procedure for dealing with it, but they 
make no mention of physical effects. Taff [Ref. 18, p 340] takes a firmer position 
by stating that it is not a prediction of second-order perturbation theory that 
there are infinitely large or infinitely rapid changes in the orbital elements, and 
that the correct resolution of these “unphysical” predictions is not to rely on bad 
mathematics. 

This present analysis lends support to TafF’s assertion. Although the criti- 
cal inclination problem was manifested in the perturbation process, it was shown 
that it does not cause singularities, nor are there any unusual physical effects pre- 
dicted by this theory. In fact, it was a remarkable aspect of this analysis that 
in the limit as the inclination approached the critical, all potential singularities 
were canceled. This is not true of the analyses which use the technique of De- 
launey/von Ziepel canonic transformations. A remarkable aspect of the canonic 
transformation method is that the non-integrability of the equations of motion 
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is not obvious when the system is expressed in coordinates or elements, but after 
the canonical transformation it is quite clear (Ref. 42] . The critical inclination is 
an intrinsic singularity of the method [Ref. 32]. 
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